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FOREWORD 


The SPAR system of computer programs was designed and developed 
by Lockheed Missiles and Space Company for the National Aeronautics and 
Space Administration under Contract NAS8-26352. This contract was funded 
jointly by NASA's Marshall and Langley Centers, and the U. S. Army 
Missile Command. The Contracting Officers ' representatives were 

L " A * ° £ the MSFC Aero-Astrodynamics Laboratory and 

J. C. Robinson of the Structures Division at Langley. 



CONTENTS 


Section 


FOREWORD 

1 INTRODUCTION 

2 GENERAL INFORMATION 

2. 1 Reference Frames. 

2.2 Card Input Rules . 

2. 3 Processor Execution. 

2.4 RESET Controls and Core Size Modification. 

2. 5 Processor Output Grouping. 

2. 6 Data Salvage. 

2. 7 Sparse Matrix Solution Method. 

2.8 Error Messages . 


3 PROCESSORS 

3. 1 TAB. Basic table input. 

3.1.1 TEXT 

3. 1. 2 MATERIAL CONSTANTS (MATC) 

3. 1. 3 DISTRIBUTED WEIGHT (NSW) 

3. 1. 4 ALTERNATE REFERENCE FRAMES (ALTREF) 
3. 1. 5 JOINT LOCATIONS (JLOC) 

3 . 1.6 JOINT REFERENCE FRAMES (JREF) 

3.1.7 BEAM ORIENTATION (MREF) 

3. 1. 8 BEAM RIGID LINKS (BRL) 

3.1.9 E21 SECTION PROPERTIES (BA) 

3.1.10 BEAM S 6x6 (BB) 

3. 1. 11 E23 SECTION PROPERTIES (BC) 

3.1.12 E24 SECTION PROPERTIES (BD) 


11 




3. 1. 13 SHELL SECTION PROPERTIES (SA) 

3. 1. 14 PANEL SECTION PROPERTIES (SB) 

3. 1. 15 CONSTRAINT DEFINITION (CON) 

3. 1. 16 JOINT ELIMINATION SEQUENCE (JSEQ) 
3.1.17 RIGID MASSES (RMASS) 

3.2 ELD. 

Element definition processor. 

3.3 TOPO. 

Element topology analyzer. 

3.4 E. 

Element information packets . 

3.5 EKS. 

Element information packets. 

3. 6 K. 

Forms system K. 

3. 7 M. 

Forms system M. 

3.8 KG. 

Forms system Kg. 

3. 9 LCM. 

Combines system matrices. 

3. 10 MRM. 

Adds rigid mass to system M,, 

3. 11 INV. 

Factors system matrices. 

3. 12 Q. 

Reads applied load data. 

3. 13 DSOL. 

Computes static solutions. 

3. 14 LCD. 

Combines solutions. 

3. 15 GSF. 

Computes stresses. 

3. 16 PSF. 

Prints stresses. 

3. 17 EIG. 

Eigensolver . 

4. DCU, THE DATA COMPLEX UTILITY PROGRAM 

5. REFERENCES 


APPENDIX A: Short List of Commands and Data Sequences. 
APPENDIX B: Element Formulation Summary. 


iii 



Section 1 
INTRODUCTION 

SPAR is a system of related programs which may be operated 
either in batch or Demand (teletype) mode. Information exchange 
between programs is automatically accomplished through one or more 
direct access libraries, known collectively as the data complex. Card 
input is command-oriented, in free-field form. Capabilities available 
in the first production release of the system include the following. 

• Highly automated data management. Data complex 
size is unlimited. 

• Linear stress analysis. Pressure, inertia, and thermal 
loading is permitted, in addition to concentrated forces 
and specified motion. 

• Linear bifurcation buckling analysis, including effects 
of prestress. 

• Linear vibrational analysis, including effects of prestress. 

Low computer execution costs and large size capacity are achieved 
using sparse matrix techniques substantially the same as those described 
in Reference (1), and implemented in several earlier programs, including' 
SNAP (Reference 2). Minimal core requirements permit the system to be 
operated economically in Demand (teletype) mode. Dynamic modification 
of core size is permitted. Most SPAR executions rarely require more 
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than 25, 000^ Q - 30, 000 l0 core locations (total of instruction + data storage). 

As an example of siee capacity, static solutions to 10,000 10 d.o.f. problems 
can be computed within a total field length of 35, 0<K> 10 . Much larger problems 

can be solved if more core is used. The size capacity of the eigensolver is 
similar. 

On UNIVAC-1100 systems SPAR is comprised of an array of separate 
absolute programs. On CDC-6000 systems, the SPAR programs appear as 
primary overlays sheltered by a minimum executive. The SPAR CDC exec 
is configured to allow the same control card/data stream to be read on both 

UNIVAC and CDC systems. That is, CDC card input records appear as 
shown below: 

@XQT PROGX 

Data cards read during execution of PROGX. 

@ XQT PROGZ 

Data cards read during execution of PROGZ. 

® XQT PROGY 
etc. 

7/8/9 

In the above example, each of the processors, PROGX, PROGZ, - - - - 
may obtain information from the data complex, and insert data into it. 

Each direct access library is a drum or disk "file" (in both UNIVAC and 

CDC terminology). Accordingly the analyst may discontinue computation 

at any point, and retain all results using either of two simple methods: 
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(1) By assigning the libraries as cataloged files using 
standard control card language provided for that 
purpose on both UNIVAC and CDC systems, or 

(2) By copying the file(s) from disk to tape using either (a) system 
commands (e. g^COPY, GM), or (b) simple copy commands 
available in SPAR's data complex utility program, DCU. 

To restart, the files are copied back to disk. 

The analyst has extensive control of which libraries are used as data 
sources for each processor, and into which libraries processor output 
dataflows. 

In problems of small-to-moderate size, the usual practice is to 
keep all data in one library (i. e. , the data complex is a single file). 
Ultimate simplicity of operation is obtained in this way, especially if the 
file is cataloged. To resume an earlier analysis, for example, it is 
sufficient to @ASG (or ATTACH) the file, then SPAR processors may be 
executed as though the previous run had not been terminated. 

Within the SPAR processors, the libraries are identified by 
logical unit numbers 1, 2, 3, - - - 26. These are not Fortran unit numbers 
The UNIVAC internal file names corresponding to libraries 1 , 2, 
are SPAR-A, SPAR-B, - - -. On CDC the corresponding file names are 

SPARLA, SPARLB, . SPAR reserves units 21-26 for temporary 

(scratch) use by processors. On both CDC and UNIVAC systems, SPAR 
dynamically assigns new files as they are needed. The analyst assigns 
(or attaches) existing fil'es, if needed. 
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Each library has imbedded in it a Table of Contents (TOC), which 
is automatically maintained by the SPAR information management system 
(IMSYS), through which processors communicate with the data complex. 
When a processor inserts a new data set into a library, IMSYS adds to the 
TOC a line containing the following information: 

• The four-word name, such as VXBR EVAL 1 0 , identifying 
data set. Usually the names are automatically assigned by 
the originating processor, although the analyst may 

occasionally assign part or all of a name through control 
card input. 

0 The date and time the data set was created. 

• Data type and error codes . If an abnormal condition (e.g. 
excessive round-off error) is detected by the processor 
originating the data set, the error code is set to indicate 
(1) warning, or (2) fatal error, as appropriate. This 
makes it possible for processors to avoid wasted effort 
by automatically declining to operate on data marked 

in error. 

• Several words describing the size and configuration of the 
data se t. A single data set may be millions of words long, 
comprised of many physical records, or it may be a small 
single -block entry, or anything in between. 

® Pis e address information enabling IMSYS to retrieve the 
data upon request. 
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TOC printouts, such as the one shown on Figure 1-i, are very 
useful in helping analysts keep track of the progress of large complicated 
studies. They comprise a complete record, in chronological order, of 
the creation of all data sets. 

Before IMS YS inserts a data set into a library, it scans the TOC and 
disables the existing data set {if present) having the same name. In TOC 
printouts, disabled data sets are indicated by a minus sign preceding "RR" f 
as shown on Figure 1-1. Disabled data sets can be re-enabled, if required, 
through the data complex utility program, DCU. 

If any of the four words of a data set name are replaced by MASK 
those words are not used by IMSYS in scanning the TOC. That is, a data 
set named GY B 7 412 could be referenced as GY MASK 7 MASK. 

Some of the primary data sets generated by SPAR processors are 
indicated below. 

• Individual tables used in constructing the basic definition 

of the structure, e.g. material constants, section properties , 
nons tructural distributed weight, joint locations, rigid 
lumped masses, element definitions in primary form, etc. 

• Detailed element information packets. Depending upon the 
type of element, each packet usually contains the following: 

(1) General information, such as connected joint 
numbers, table reference pointers, etc. 
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(2) Material constants. 

(3) Detailed geometrical data, e. g 0 dimensions, 
transformations . 

(4) Section properties. 

(5) Intrinsic stiffness matrix, or equivalent. 

(6) Intrinsic stress matrix, or equivalent. 

® System M, K, Kg, matrices, and linear combinations of 
them, in SPAR's sparse matrix format. 

• Factored system matrices in sparse matrix format. 

• Applied loading and constraint data. 

• Solution data such as displacements, stresses, 
eigensolutions , etc. 
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Section 2 


GENERAL. INFORMATION 

It is assumed that the reader is familiar with linear finite element 
theory, as applied to static, vibrational, and bifurcation buckling analysis 
of structures. 


Table ELD-1, near the beginning of Section 3.2, contains a brief 
summary of the types of elements currently available in SPAR. Appendix 
B gives details of the formulations on which individual element K, M, and 
Kg calculations are based. Reference 3 contains a series of annotated 
example problems. 

Before reading further, the new user should examine Table ELD-1, 
and briefly scan Reference 3 to get a general idea of SPAR’s control 
card - data card structure. 


2-1 



2. 1 Reference Frames 


The term frame k will be used to refer to the global reference 
frame <k = 1), or to any alternate reference frame (k = 2, 3, 4, - - ) the 
analyst elects to define. Each joint in the structure has associated with 
it a unique joint reference frame , to which joint displacement and rotation 
components are relative. Individual joint reference frames may have any 
orientation designated by the analyst. Each element has an associated 
element refere nce frame , to which section properties and stress 
components are relative. 
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2. 2 Card Input Rules 

The SPAR free -form input decoder recognizes three types of 
words: integer, floating-point, and typeless. Leading blanks are 
ignored. Each word is ended by a blank, comma , equal sign , slash, 
or by a record terminator (e.g. end-of-card). If used, commas, etc. , 
should be carefully placed; for example 45 , is equivalent to 45, 0, . 
Floating-point numbers are identified by the presence of a decimal 
point. Typeless words must begin with a letter. Allowable forms for 
floating-point numbers are: 

5000.= . 5+4= . 05+5= 50.00+02, etc. 

Note that the Fortran form x.xExx is not permitted. 

ft 

Each card begins an input record. A record is terminated by 
end-of-card, or by any of the three following symbols: 

Record 

T erminator Function 

$ Characters to the right of a $ are ignored. This is 

used for two purposes: (1) to allow the input decoder to 
stop scanning, and (2) to allow the user to insert comments 
in the data deck. A card with a $ in Column 1 is 
interpreted as a comment card, and is ignored by the 
decoder. 

5/8 punch The 5/8 punch ( : on Univac, < on CDC ) terminates one 
record and initiates a new record on the same card. 

For example, 


2-3 



3, 4 : A, 9.5 $ 
is the same as 

3. 4 $ 

A, 9.5 $ 

4/8 Punch All characters to the right of a 4/8 punch { > on Univac, 

< on CDC) form a continuous alphanumeric label of up 

to 76 characters . Examples of label usage are shown 
below. 

GROUP 42' RING FRAME , STA 0 420. 

CASE 240' 2, 5 G GUST + 9.4 PSI PRESSURE. 

Throughout this document, in put card descriptions will use frh* 
Umvac symbols ( : for 5/8 pun ch, ' for 4/8 punch, ® for 7/8 punch) . 

Typeless words longer than 4 characters are truncated to 4 
characters. Integer or floating-point words must not exceed 7 digits. 

Exponents may have one or two digits. Floating-point numbers must not 
exceed host system limits. 
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2.3 


Processor Execution 


The short summary of processor functions given at the beginning 
of Section 3 should be thoroughly examined before reading further in the 
present discussion. 


Typical execution sequences are illustrated below. 


Basic 

Preparations 

@ XQT TAB 
input cards 

@ XQT ELD 
input cards 

@ XQT TOPO 

@ XQT E 

@ XQT EKS 

® XQT K 

@ XQT M; if 
required 

@ XQT INV 


Static Buckling or 

Solutions Vibrational Modes 


@ XQT Q @ XQT EIG 

input cards control cards 

@ XQT DSOL 


@ XQT LCD 
input cards 

@ XQT GSF 

@ XQT PSF 


Form System Kg 
@ XQT KG 


Of the above, only TAB, ELD, and Q read significant amounts 
of data through input cards. In TAB, the analyst is permitted to make 
corrections to individual items in tables of joint locations, section and 
material properties, etc. This is done by modifying previously created 
data sets resident in direct access libraries, thereby greatly reducing 
the card input load on the host system. It is expected that similar pin- 
point editing capability will be released in the near future for data sets 
produced by ELD, Q, and other processors. 


2-5 



In preparing finite element models of complicated structures, 
it is almost invariably best to begin by preparing separate models of 
individual components and analyzing them to see that they perform 
reasonably; e.g. , by computing response to representative static 
loadings. Models of components are then merged, perhaps in several 
stages, with additional testing performed after each merger. In this 
way, input errors are easily and systematically identified. The analyst 
working on each individual model usually finds it convenient to use 
reference frames, joint numbers, section property tables, etc. , local to 
the particular component. In TAB and ELD, various data modification 
commands are provided to assist in merging models. It is suggested 
that the new user carefully examine the NREF and MOD commands in 
TAB, and the MOD and INC commands in ELD. Additional devices to 
further automate the merging of models are planned for future releases. 

Most of the processors generate little or no printout. In the few 
that do, some degree of control of the quantity is provided by inserting 
a data card containing the following command: 

ONLINE = n, 

where n = 0 for minimum (or no ) printout, n = 1 for normal printout, 
n = 2 for more. 
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Reset Controls and Core Size Modification 


RESET statements allow the analyst to alter certain parameters 
which control processor operation. Typical parameters available through 
RESET are listed below. 

Names of input and output data sets. 

Unit numbers of source and destination libraries 
(default is always SPAR logical 1). 

Flags indicating action to be taken if errors are encountered. 
Flags controlling display models . 

Miscellaneous quantities, such as zero-test values, scale 
factors, output data set block lengths , formulation selections , 
constants, iteration controls, etc. 

The form of a RESET statement is: 

RESET P 2 =w 2> p 3 =w 3> , where w. is the value to be 

assigned to parameter p. . One or more RESET statements may immediately 
follow a ®XQT P rocessor command : for example: 

@ XQT INV 

RESET K=K+KG, CON =2 $ 

® XQT EIG 

RESET PROBLEM = STABILITY! 

RESET CON=2, INIT=7, NDYN=8 $ 

In addition to specific RESET controls defined for individual 
processors, the following two parameters apply to all processors: 

RESET ABORT = 1 
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which will cause the processor not to make an error abort if it encounters 
a serious error (e.g. , if required input data sets do not exist), 

RESET CORE = n (available on UNIVAC, only), 

which will result in issuance of an executive request to change the total 
core size (both instructions and data) to n words. 

On CDC systems, the user controls core size through RFL cards. 


2-8 



Processor Output Grouping 


The RESET controls and the data complex utility program, DCU, 
provide the analyst total control over the distribution of data sets among 
libraries. However, the following rules should be observed uniformly. 

• Output from TAB, ELD, E, and EKS should always be in 
Library 1. Library 1 should be assigned when any processor 
(except DCU) is executed, since it contains several small 
directories which are widely used. 

• For each associated set of quasi-static load cases, the 
output of Q, DSQ L, GSF and LCD (if us ed) must be in the 
same library. It is usually best not to put more than a 
few dozen independent loadcases into a single library. 

Since GSF output data sets are comparatively large, it is 
often best to save only the Q and DSOL data sets on archive 
tapes, in order to cut down on the amount of stored data. 
Stresses can be recalculated very cheaply, if necessary. 
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2. 6 Data Salvage 

If it is necessary to make limited changes in the model (e.g. 
to study effects of parameter variations, or to patch errors), it is 
often possible to reuse much of the existing data complex, provided the 
analyst is sufficiently aware of how information flows among the processors 
(see Figure 2-1). A few examples are illustrated below: 

• Vibrational modes have been computed for a given 
configuration. It is now required to compute modes 
for the same structure with different rigid masses 
attached. The procedure is: (1) re-execute TAB sub- 
processor RMASS to read the new rigid mass data, (2) 
re-execute MRM to add the new mass data to the 
existing system M representing element mass, and (3) 
re-execute EIG, using the previously computed modes 
as initial approximations. 

• After a system matrix has been factored, it is discovered 
that incorrect constraint conditions were used. The 
procedure is: (1) re-execute TAB sub-processor CON to 
patch the constraint case definition, and (2) re-execute 
INV. 

• After computing vibrational modes, it is discovered that 
there is an error in cross -section property data resulting 
m a bad system K, but not affecting M. The procedure is 
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(1) re-execute the appropriate TAB sub-processor to 
correct the erroneous item, (2) re-execute E, EKS, 

K, IKV, and (3) re-execute EIQ, using the previously 

computed modes as initial eigenvector approximations, 
if appropriate. 


In selecting processor execution sequences in cases such as 


those indicated above, users who 
should proceed with caution. 


are not familiar with the system 
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2 * 7 Sparse Matrix Solution Method 

Reference 1 presents a procedure for solving high-order systems 
of linear equations of the kind which occur in displacement-method finite 
element analysis. The system stiffness matrix is regarded as an array 
of submatrices, each submatrix corresponding to the connection of one 
joint to another. Each submatrix is n by n, where n is the number of 
degrees of freedom at each joint. For general shells and frames, n = 6 
(three displacements, three rotations); for a plane frame, n = 3 (two 
displacements, one rotation); etc. The only nonzero submatrices in a 
system stiffness matrix are those corresponding to pairs of joints con- 
nected by elements. Accordingly in all but the smallest finite element 
models only a tiny fraction of the submatrices are nonzero - usually less 
than 1%. The characterizing feature of the procedure described in 
Reference 1 is that it operates exclusively with data contained in the non- 
zero submatrices, virtually eliminating the unessential arithmetic 
(multiplying, adding zeros) and wasted data storage (storing zeros) 
associated with conventional band matrix tec hni ques 

Originated in the early 1960's this method was first applied in a 
series of large capacity frame analyzers, and later in the SNAP series 
of finite element network analyzers. SPAR uses essentially the same 
procedure, except for certain differences in the order in which certain of 
its phases are performed. It will be necessary for users of SPAR to become 
familiar with Reference 1 ; especially with regard to the effects of joint 
numbering sequence on execution costs and data storage requirements. 
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2. 8 Error Messages 

Various error messages, usually self-explanatory, may originate 
in SPAR processors. Messages indicating insufficiait core space refer 
to working core requirements, not to total (instruction + all data) space. 

At the beginning of execution of each processor, a message, DATA SPACE = 
XXXXX, indicates the amount of working core space currently available. 

This information enables the analyst to estimate the total field length 
needed in the rerun. 

Another common message is ERROR IN CONTROL, STATEMENT 
SYNTAX. This means that an illegal RESET name has been given, or that 
Ihe wrong type of data (integer, floating-point, typeless) was given. 

The most frequently encountered error messages indicate that 
either a required input data set is not present (IERR=-1), or that it is 

marked m error (IERR = a positive integr). Both the causes and remedies 
for such messages are usually obvious. 


The messages XL, IB - n and ASG = n indicate an attempt to refer 
to a nonexistent library, n. 
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Section 3 


PROCESSORS 


SPAR processors currently available are summarized below. 
Processor 

Section Name Function 


3.1 TAB 


3. 2 ELD 


Through, an array of sub-processors, card 
input is translated into tables of material 
constants, section properties, joint locations, 
rigid mass data, constraint conditions (any 
number of separate cases are allowed), 
joint elimination sequence, etc. 

Element definitions are read from cards. 

The definitions refer to entries in the TAB- 
generated tables of section and material 
constants, etc. Various automatic mesh 
generation options are provided. 


Element interconnection is analyzed and two 
data sets produced: one guides assembly of 
system K, M, and Kg matrices in SPAR's 
sparse matrix format; the other guides factoring 
of sparse matrices by processor INV. Size and 
cost indices are computed permitting accurate 
estimates to be made of the execution costs which 
will be required by processors INV and DSOL. 
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Processors (continued) 


Section 

3.4 


3 . 5 


3 . 6 
3 . 7 

3.8 

3.9 


3 . 10 


Processor 

Name Function 

E Using data sets created by TAB and EL»D, 

information packets are formed to describe 
each element in detail. The packets contain 
geometrical details, section and material 
properties, and other information. 


EKS 


K 

M 

KG 

LCM 


MRM 


Element intrinsic stiffness and stress matrices 
(or their equivalent) are computed and inserted 
into the element information packets. These 
matrices are based on current values of section 
properties, etc. , imbedded in the element 
information packets. 

System K is formed in sparse matrix format. 
System M is formed in sparse matrix format. 
System Kg is formed in sparse matrix format. 
Any required linear combinations of system 
matrices in sparse format are formed. 
Examples: K+Kg, K-rM, aKj+bK^ etc. 

Adds a diagonal matrix to a SPAR-format 
matrix. Normally used to add effects of rigid 
mass to a system M. 
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Processors (continued) 


Processor 
Section Name 

3.11 INV 


3. 12 Q 


3.13 DSOL 


Function 

SPAR-format system matrices such as K, 

K+Kg, K-rM, etc. , are factored using a 
sparse matrix procedure substantially the 
same as the one described in Reference 1. 

Through sub-processors, card input is read 
defining applied quasi- static loading. The 
following types of loading are included: 

• Inertia loads 

• Pressure distribution 

• Temperature distribution 

• Specified nodal motions 

• Concentrated loads. 

Q creates data sets defining load case details, 
but does not execute the corresponding solutions. 

Solutions (joint motions) due to selected load 
cases are computed. 


3.14 LCD 


Linear combinations of (joint motion) 
solutions are formed. 


3.15 GSF 


Data sets containing element stresses and 
internal loads are generated. 


3. 16 PSF 


Stress data printout is generated. 
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Processors (continued) 


Processor 
Section Name 


Function 


Vibrational and buckling eigenproblems are 
solved. High-order problems represented 
by SPAR format M, K, Kg, etc. , matrices 
are solved directly. 


The data complex utility program, used to 
copy data sets from one library to another, 
change data set names , disable ( 0 r re-enable) 
data sets, print data set contents, print 
library tables of contents, etc. 
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3. 1 TAB 

s 

Function. TAB contains an array of sub-processors which 
generate tables of material constants , section properties, joint 
locations, etc., and. various other data sets comprising a substantial 
portion of the definition of the structure. Each of these sub-processors 
is identified either by a multi-word (long form) name, such as JOINT 
LOCATIONS, or by a short name, e. g, , JLOC. Each sub-processor 
generates a data set having the same name as the sub-processor. 
Sub-processor names and functions are summarized in Table TAB-1. 

TAB may be used to either (1) create new data sets, or (2) 
update existing data sets by replacing individual entries in them. The 
update mode is commonly used in Demand (teletype) interactive operation. 
Only one direct access library, SPAR logical 1, is usable in conjunction 
with TAB. When beginning a new problem , the first data card following 
@ XQT TAB must be 

START j, m^, m^, - - - 

In the above, j = the total number of joints in the structure. It is not 
harmful to have some unused joints (i. e. joints connected to no elements), 
for convenience in interpreting the output. This should not be carried to 
extremes, however, since it wastes core storage. It should be noted that 
it is not necessary to set aside large blocks of unused joint numbers in 
areas where you expect to later augment the model. The JSEQ sub- 
processor and an array of data modifier statements allow models to be 

extended without extensive repunching of existing data cards. 
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The parameters on the START card list joint 

motion components which are identically zero for all joints. Those 
components are relative to the joint reference frames. For m=l,2, or 3 a 
direction m displacement is indicated. For m=4, 5, or 6, a direction 
m-3 rotation is indicated. 

For example, a 1000 joint space truss would begin with 
START 1000, 4,5,6. 

The above indicates that all three rotation components are zero at all 
joints. As another example, a 500 joint plane frame would begin with 

START 500 3, 4, 5. 

The above assumes that the the direction-3 displacements are normal 
to the plane of the frame, and the frame deforms in-plane only. 

SPAR uses the m's given on the START card to determine 
submatrix storage block size (see Reference 1). Therefore the START 
card (rather than the CON sub-processor) should be used to suppress 
joint motions in cases such as those illustrated above. 

Following the START card, or directly following @ XQT TAB if 

Library 1 already exists, the TAB card input stream is as illustrated 
below. 

Proc^ 
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data cards read by Proc 


1 


Proc 2 

data cards read by Proc 

2 


etc. 

In the above, Proc. represents a data card containing either the short 
or long-form name of a sub-processor to be executed. Sub -processors 
may be executed in any order, subject to the following restriction: if 
the data read by one sub-processor refers to another TAB-generated 
table, the other table must already exist. It is always safe to execute 
sub-processors in the order they appear in Table TAB-1. 

TAB input card rules 

* M addition to the rules generally applying to SPAR free -field 
input, the following rules apply specifically to TAB. 

(1) Trailing items omitted on data cards are assumed to be 
zero, except when indicated otherwise in discussions of 
specific sub-processors • 

(2) Both integer and real data may be given in either fixed 

or floating point format (i. e. the integer 5 may be given as 5. , 
and the real variable 5. may be given as 5 ). 
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(3) Within, the data card stream being read by any sub— processor, 
the following commands may be injected. 

FORMAT = j 
MOD= m 

NREF= n 

In the above, 

• j identifies which of several admissible 
formats is to apply to subsequent cards. 

Details are given in discussions of individual 
sub-processors. 

• m is added to the table entry indicator 
(usually called k), in all processors other than 
TEXT, JLOC, JREF, CON and JSEQ. In JLOC, 

JREF , and CON, m is added to joint numbers 
indicated on subsequent input cards. 

• Ln processors JLOC, JREF, and MREF, reference 
frame n applies to data on subsequent cards. 


These three commands may be used repeatedly in the card 
streamof any sub-processor. Upon beginning execution of a 
new sub-processor, these parameters are internally reset to 
their default values, FORMAT=l, MOD= 0 , and NREF=1. 


Correcting existing data sets 


To enter the "update" mode of operation, the following command is 
injected in the input stream. 


3. 1-4 



TAB 


UPDATE = 1 

To leave the update mode, the command is 
UPDATE =0 

UPDATE commands should immediately precede: sub -processor execution 
commands. When operating in the update mode, the output data set produced 
in the current execution is identical to that produced in the preceding 
execution, except for entries defined by the card input of the current 
execution. The only sub-processors which cannot be operated in the 
update mode are TEXT and JSEQ. 

As one example, suppose the location of joint 1742 is found 
to be in error. The JLOC data set could be repaired by the following 
card sequence. 

UPDATE=1 

JLOC 

1742, 947.62, 1841.9 23.487$ 

Reset Controls 

No special RESET controls are available in TAB. 

Core Requirements 

Working core requirements will not significantly exceed the 
larger of (1) 13 times the number of joints, or (2) the longest table 
generated (see Table TAB-1). 

Code Release Data 

Level 7 (UNI VAC, CDC) July 1974, coded by W. D. Whetstone, 
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Table TAB-1, TAB Sub-Processors 


Short 

Name 

Table . 
Length 

TEXT 

20 n 

matc 

10 n 

NSW 

n 

ALTREF 

12 n 


Long Form Name, 

Function 

TEXT. Creates a data set containing 
alphanumeric text documenting the analysis. 

MATERIAL CONSTANTS. Entries define 
linear material constants E, p, etc. 

DISTRIBUTED WEIGHT. Entries define 
nonstructural distributed weight parameters 
(weight/ length or weight/area). 

ALTERNATE REFERENCE FRAMES. Entries 
describe reference frames selected by the 

analyst for convenience in defining joint locations, 
etc. 


JLOC 


3 j JOINT LOCATIONS. 


j JOINT REFERENCE FRAMES. Entries define 

the orientation of reference frames associated 
with joints, used in defining constraint, applied 
loading, etc. 


MREF 

5 n 

BRL 

8 n 


BA 

31 n 

BB 

21 n 


BEAM ORIENTATION. Entries define beam . ■ 
cross-section orientation. 

BEAM RIGID LINKS. Entries define rigid 
links offsetting the end points of elastic 
2 -node elements from the joints they connect. 

E21 SECTION PROPERTIES. 

BEAM S 6 x 6 . Entries define elastic 
characteristics of type E22 and E25 elements. 


* n is the number o f entries in the table, and j is the total number 
of joints in the structure. 
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Table TAB- 

-1, continued 

BC 

6 n 

E23 SECTION PROPERTIES. 

BD 

6 n 

E24 SECTION PROPERTIES. 

SA 

19 n 

SHELL SECTION PROPERTIES. 

SB 

4 n 

PANEL SECTION PROPERTIES. 

CON 

j, 6j 

CONSTRAINT DEFINITION. 

JSEQ 

j 

JOINT ELIMINATION SEQUENCE 

RMASS 

6j 

RIGID MASSES. 
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3.1.1 TEXT 


TAB 
TEXT 
MATC 
NSW 

The function of the TEXT sub-processor is to allow the analyst 
to imbed in the output library a data set containing alphanumeric text 
descriptive of the analysis being performed. Each card has a 4/8 punch 
in column one, followed by a 76 character string of text. 

3 * 1 * 2 MATERIAL CONSTANTS (MATC) 

MATC generates a table of linear material constants. Each 
input card defines one entry in the table. The data sequence on the input 

card defining the k-th entry in the table is k, E, Nu, Rho, Alpha. 

In the above, 

E = modulus of elasticity, 

Nu = Poisson's ratio, 

Rho = Weight per unit volume, and 
Alpha = Linear thermal expansion coefficient. 

Material constants associated with specific elements are defined 
in processor ELD by pointing to entries in the MATC table. 

3.1.3 DISTRIBUTED WEIGHT (NSW) 

A table of non- structural distributed weight parameters is defined. 
The data sequence of the input card defining the k-th entry in the table is 
k, W. 

For 2-node elements, W is weight/length. 

For 3 and 4 node elements, W is weight/area. 

Non- structural weight attached to specific elements is defined in 
processor ELD, by pointing to entries in the NSW table. 


3. 1-8 



tab 

ALTREF 


3 .i.4 ALTERNATE REFERENCE FRAMES (ALTREF) 

In addition to the global reference frame, the analyst may find 
it convenient to define additional reference frames. These frames have 
several uses, including the following: 

(1) Joint locations may be defined in any frame the analyst 
finds most convenient (see JLOC). 

(2) Joint reference frame orientations may be defined via 
the alternate frames (see JREF). 

•^ ac h frame is uniquely identified by a positive integer. The 
global fr ame i s always frame 1 ; accordingly the analyst is free to define 

only frames 2, 3, The order of data on the input card defining 

frame k is 

k, i 1 ,a r i 2 ,a 2 , i 3 ,a 3 , x v x 2 , x^ 

In the above Xj, and x 3 are position coordinates, relative to the 
global frame, of the origin of frame k. The x's need not be given if 
only the orientation of frame k is of significance, which often is the case. 

The parameters i Jf a 3 indicate ordered rotations 

defining axis orientations. Two formats are provided. In both cases, we 
begin with a local frame parallel to the global frame. After the ordered 
rotations are completed, the local frame is parallel to frame k. Each of 
the i's may be 1, 2, or 3, in any order. The a's are angles in degrees . 

If FORMAT=l (default) , the sequence is: (1) rotate the local 
frame &1 degrees about local axis then (2) from the new position. 


3. 1.4-1 



TAB 

ALTREF 


rotate the local frame a 2 degrees about axis i^, then (3) from the resulting 
position, rotate the local frame a 3 degrees about axis i y If FORMAT =2, 
the i's and a's indicate rotation of the global frame relative to frame k. 

In the following example of an input card defining frame 27, it 
is assumed that FORMAT = l. 

27 3,30. 1,10.$ 



G. - Global axes 

1 K. - axes of Frame 27. 

i 
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3.1.5 JOINT LOCATIONS (JLOC) 

JLOC produces a table containing the position coordinates of 
the joints. The data sequence on input cards is as follows: 

k, xa^, xa^, xa.^, xbj, xb 2 , xb^, ni, ijump, nj 

If nj is given, a second card must appear, 

jjump, xc ^ , xc ^ » xc^, xdj, xd 2 , xd^ 

There are three possible interpretations of the above: 

(1) If only k, xaj, xa^ and xa^ are given, the xa's are coordinates 
of joint k. 


( 2 ) 


If k, xa^, - - - - xbj, ni, and ijump are given, the xa's and 
xb's are coordinates of points A and B terminating a string of 
ni equally-spaced joints, 


k, k+ijump, k+2(ijump) - - 


Point B 
* 


Point 


A 

[ 

k 


k+ 

ijump 

The default value of ijump is 


k+ 

2 (ijump) 


1. 


k+ 

(ni-1) (ijump) 


(3) If nj is given, a linearly interpolated two-dimensional mesh of 
(ni)(nj) joints are defined, as indicated below. 


C 

+ 

k+ 

2 (jjump) -+- + 

k+ 

jjump 


A + 
k 


+ 

+ 




■Jr 




+ 


-V 


A- 


k+ 

-f- ijump+ + 
jjump 


+ 

k+ 

ijump 


-t- 

k+ 

2 (ijump) 


4 


At 


-VB 
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Although the output table generated by JL.OC is in rectangular 
coordinates relative to the global frame, coordinate data appearing on 
the input cards may be either rectangular or cylindrical, and may be 
relative to any frame already defined in ALTREF. The associated 
control cards are summarized below. 

Control Meaning 

NREF = n Coordinate data on subsequent cards is relative 
to frame n (until another NREF command is 
encountered). 


FORMAT- 2 Subsequent data is in cylindrical coordinates, 
relative either to frame 1 (global) or to any 
other frame selected by an NREF command. 
The convention is shown below. 



X 1 ~ ra dial distance from 3 -axis. 

x t - angle, in degrees , from 
the 1 - 3 plane . 

x 3 = linear distance from the 
1-2 plane . 


FORMA T - 1 Switch back to rectangular coordinates. 


Switching among frames, and between rectangular and cylindrical 
coordinates is unrestricted. 


If cylmdrical coordinates are used in connection with mesh generation, 
interpolation is performed before transformation to rectangular coordinates; 
so that regular meshes on circles, cylinder and cones are readily generated. 
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3.1.6 JOINT REFERENCE FRAMES (JREF) 

A unique right-hand rectangular reference frame is associated 
with each joint. Through JREF the analyst may designate the orientation 
of any joint reference frames. All joint reference frames not defined in 
JREF are, by default, parallel to frame 1 (global). The orientation of 
these frames is of considerable significance, since joint motion components 
and mechanical loads applied at joints are defined relative to them (see 
discussions of TAB sub-processor CON, and processor Q). 

To identify a set of joints with reference frames parallel to 
frame n (see ALTREF), an input card containing NREF=n is followed by 
a sequence of records in 'loop limit 1 format naming the joints. 'Loop 
limit' format is: jj, j 2> inc, meaning joints j^inc, jj+2(inc), j^. 

Default for j 2 is j . Default for inc is 1. 

For example, the following records indicate that the reference 
frames of joints 5, 6, 7, 8, 19, 30, 32, 34, and 36 are parallel to frame 10. 
NREF=10: 5,8: 19 : 30,36,2$ 

The language NREF=-n indicates that the joint reference frame is 
oriented as follows: 

• The 3-axis of the joint frame is parallel to the 3-axis of 
frame n 
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• The 1 axis of the joint frame is perpendicular to the 
3-axis of frame-n, as shown below : 
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3*1*7 BEAM ORIENTATION (MREF) 

Each two-node element has an "element reference frame" 
associated with it (see processor ELD discussion). Beam section 
properties, stresses, etc., are defined relative to these frames. 

The 3 -axis of the frame is directed from the beam’s origin to its 
terminus. The origin and terminus are the ends connected to joints 
J1 and J2, respectively, as defined in processor ELD. The beams 

end points coincide with J1 and J2, unless offset by rigid links defined 
via the BRL sub-processor . 

Entries in the table produced by MREF are used to define 
the orientation of beam axis 1 and 2. Usually a single entry in this 
table will apply to many beams. Either of two data sequences may be 
selected through FORMAT control, as indicated below. 

r 

If F ORMA T=1 (default), the data sequence defining table entry 

k is 

k, nb, ng, isign, c 

The above indicates that c is the cosine of the (smallest) angle between 
beam axis nb and global axis ng. Legal values of nb are 1 and 2, and 
ng may be l,2,or3. The parameter isign resolves a possible ambiguity 
by indicating whether the cosine of the angle between beam axis (3-nb) and 
global axis ng is positive or negative, with values of +1 and -1 indicating 
positive and negative, respectively. Examples are shown below. 
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nb, ng, isign, and c are nb, ng, isign, and c are 

1, 3, -1, . 1. 1, 3, +1, . 1. 

In both of the above cases the cosine of the angle between beam axis 1 

and global axis 3 is .1, and isign defines the sign of the cosine of the 

angle between beam axis 2 and global axis 3. 


TAB 

MREF 


In about 95% of the cases c is either 1.0 or . 0 o 


Care should be taken to supply meaningful information. For 
example, if the beam 3 axis is parallel to global axis 2, it is useless 
to state that the cosine of the angle between global axis 2 and beam axis 
1 (or 2 ) is zero. 


In FORMAT-2 , the data sequence defining table entry K is 
k, il, X| , x^, x^ . 

The above indicates that beam axis orientation is determined by the 
location of a 'third point 1 , located at (x^, x^, x^)> as illustrated below. 
The x's are in rectangular components, relative to any frame defined 
in ALTREF. The NREF=n command may be used to switch to frame n 
(default is frame 1), anywhere in the input stream. 
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Origin 

In the above, E^, E 3 and B are coplanar, and E 2 is parallel to E^ x B. . 
The input parameter il selects one of four possible beam axis orientations, 
as shown below. 


E2 E2 Ej E2 



The analyst will find both formats 1 and 2 useful. For example, 
if rings and stringers reinforce a cylindrical or conical shell, it may be 
possible to define the orientation of all of the ring elements with a single 
format- 1 MREF table entry, and to define the orientation of all of the 
stringer elements with one format-2 entry. 
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3.1.8 BEAM RIGID LINKS (BRL) 

Each entry in the table generated by BRL defines rigid link offsets 
for both the origin and terminus of a 2-node element. The data sequence 
on the input card defining the k-th entry in the table is 



The above data is defined as follows. 

1 1 1 

# x^ , x^ are the position coordinates of joint 
J1 (see ELD discussion) minus the position 
coordinates of the beam origin. 


X l* x 2’ x 3 are the position coordinates of joint J 2 
minus the position coordinates of the beam terminus. 


The x's are always given in rectangular coordinates relative to reference 
frames defined by the parameters i 1 and i 2 . If i^n (positive), it indicates 

that Xj , x 2> x* are relative to frame n (see ALTREF). If i^-n the 

reference frame of x^, x 2 , x^ is shown below. 
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EZ1 SECTION PROPERTIES (BA) 

BA generates a table of section properties to which reference is 
made during definition of type E21 elements in processor ELD. Nine 
classes of cross section are allowed. 

The first word in each input record is a typeless word {e.g. BOX, 
TEE, etc. ) identifying the class. A single card defines a table entry, 

except for DSY which re quires two cards . Data sequences are as indicated 
below. 


BOX 

k. 

b i ,t r h z ,t z 

TEE 

k, 

W b 2’ fc 2 

ANG 

k, 

b l’V b 2’ t 2 

WFL 

k. 

b l ,t l’ b 2’ t 2’ b 3 ,t: 3 

CHN 

k, 

b l ,t l* h 2 f t 2’ b 3» 4 3 

ZEE 

k, 

b l * fc l ’ b 2’ t 2* b 3» 4 3 

TUBE 

k, 

inner radius, outer radius 

GIVN 

k, 

h 9 *! 9 ^2* ** 2* a * 2 1 * 

DSY 

k, 



V 

q 2* q 3* y ll’ y 12' ■ ■ “ y 4i 


(card 1) 
1 y a'? / 


In the above, k identifies the table entry number. The b's and t's are 
cross-section dimensions defined on Figure BA-1. 


In all cases the origin and terminus of the beam (see discussions 
of MREF, BRL, and ELD) coincide with the section centroid. For GIVN 
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and DSY sections, 


Principal moments of inertia. For DSY sections, principal 
axes must coincide with the element reference frame axes. 


a. , af^ Transverse shear deflection constants associated with 1^ 
and I_, respectively. For no shear deflection, set ot 
equaFto zero. 

* 

f = Uniform torsion constant. For uniform torsion, torque = 

Gf x (twist angle /unit length), where G is the shear modulus. 

fj= Nonuniform torsion constant, 0 accounting for flange -bending 

effects on torsional stiffness, etc. 

Zj,z£ Shear center - centroid offsets. 

9 Inclination of principal axes relative to the element 

reference frame (see Figure BA-1). 


Use of the above quantities in calculating element stiffness matrices 
is discussed in Appendix B. Items given on the second card defining DSY 
sections have the following meaning: 

qi.q 2 Section shape factor such that maximum shear stress due 
to V a shear in direction 1, is = V q^. q^ is similarly 
defined. 

q^ Section shape factor such that maximum stress due to 

twisting moment T is Tq^. 

^il , ^i2 = Location, relative to the element reference frame, of the 

^ i-th point at which My/I combined bending stresses are to 

be computed. Up to 4 such points may be prescribed. 


In the notation of Appendix B, 
C = Gf, and 



Ef 


1 
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3.1.10 BEAM S 6 x 6 (BB) 

BB generates a table of directly specified 6 by 6 intrinsic stiffness 
matrices to which reference is made during definition of type E22 and 
E25 elements in processor EED. Six cards , as indicated below, are 
required to define one entry in the table. 


S 11 




S 21 ’ 

S 22 



S 3 1 ’ 

s 32’ 

s 33 


S 4V 

S 42 

S 43’ 

S 44 

S 5V 

S 52’ 

S 53’ 

S 54’ 

s 6r 

s 62’ 

s 63’ 

S 64’ 


In the above, k identifies the table entry, and the , are elements 
of a symmetric intrinsic stiffness matrix, S , which is defined as follows: 

ifj he terminus (see ELD, M R EF, BRL discussions) of the element i. 
completely fixed. 


where 


F= SU, 

F = (fj f 2 f 3 rrij m 2 m 3 ) t , and 

U = (u l u 2 u 3 r l r 2 r 3 ) t * 


In the above, f., m., u., r. = appiied force, applied moment, displacement, 
and rotation components at the origin. All components are relative to the 
element reference frame. 
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3.1.11 E23 SECTION PROPERTIES (BC) 


generates a table of axial element section constants to which 
reference is made during definition of type E23 elements in processor 
ELD. The data sequence on the input card defining the k-th entry in the 
table is 


k, 


V 



max 


a . . 
mm 


In the above, 



s = 
c 


a 

max 

a . 
min 


area effective if in tension, 

area effective if in compression, 

allowable tensile stress, 

allowable compressive stress, 

maximum area allowed during redesign, and 

minimum area allowed during redesign. 


In SPAR level 7, no redesign or iterative analysis processors were 
released. Element cross-sectional area=a t . The other parameters 
a c > S £> etc. are not used. 
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3 « 1 > 12 E24 SECTION PROPERTIES (BD) 

BD generates a table of plane beam element section properties 
to which reference is made during definition of E24 elements in 
processor ELD. The data sequence on the input card defining the 
k-th entry in the table is 

k, A, Ij, <*j, h, c, q . 

In the above, 1 ^ orj and q 2 are the same as for DSY sections, as defined 

m BA, the E21 section property sub-processor. Dimensions h and c 

are as shown below. A = cross sectional area. 


Axis 2 



Axis 1 of element 
reference frame 
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3. i.13 SHELL SECTION PROPERTIES (SA) 

SA generates a table of shell section properties to which reference 
is made during definition element types E31, E32, E33, E41, E42, and 
E43, in processor ELD. Two formats are available to define table entries. 
Material constant table (MATC) entries associated with individual elements 
are defined in ELD via the NMAT parameter. Applicable entries in the 
material constant table must also be identified during execution of SA, 
using the same type of command, NMAT=en±ry number (defaults). 

Two formats are available to define entries in the SA table. If 
FORMAT=l (default), an isotropic section is indicated, and the data 
sequence on the single input card defining the k-th entry in the table is 

k> t ’ t ll 5 *22’ x iZ* 

In the above, t=thickness, and the constants t n> t 22 , and t J2 have the 
same meaning as defined below for aeolotropic sections. 

If FORMAT=2, an aeolotropic section is indicated. Three cards, 

containing the data sequences indicated below, are required to define the 
k-th entry. 

k> hi' *22’ t 12 

C ir c 21’ c 22* c 31' c 32' C 33 
d ir d 2r d 22’ d 31’ d 32' d 33 * 
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SA 




Where s j^> S 2Z‘ an< * S 12 are membrane stresses, 

S 1 1 " N ll y,t ll’ s 22~ N 22^ t 22’ S l2 =N 12^ t l2* 

For isotropic sections, default t ii =t 22 =t l2 =t ' 

AH of the above quantities are relative to the element reference 
frame (see ELD discussion). 
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3. 1. 14 PANEL. SECTION PROPERTIES (SB) 


SB generates a table of shear panel section properties to which, 
reference is made during definition of E44 elements in processor ELD. 
The data sequence on the input card defining the k-th entry is 


k. t. 


s, 


max 


min. 


In the above, 

t= panel thickness, initial design. 
s= allowable shear stress. 

t max = maximum thickness permitted during redesign. 
*min = minimum thickness permitted during redesign. 


In SPAR level 7 no redesign processors were released, and the 

parameters s, t , and t are not used, 
max mm 
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GON 

3. 1.15 CONSTRAINT DEFINITION (CON) 

Any number of constraint cases may be defined, each uniquely . 
identified by an integer assigned by the analyst. To define constraint 

case n, the CON sub-processor may be activated by any of the following 
commands . 

CONSTRAINT DEFINITION n 

CONSTRAINT = n 

CON=n 

Constraint is defined by specifying certain joint motion components 
to be either {1) unconditionally equal to zero, or (2) or specified value, 
possibly not zero. In the second category, specific values of such 

components may vary from one load case to another, as defined in 
processor Q. 

In the following, m=l, 2, or 3 indicates joint displacement in 
direction m , relative to a joint reference frame (see JREF). If 
m=4, 5, or 6, a direction m-3 rotation component is indicated. 

To specify that joint motion components m., m ? , - - are 
unconditionally zero at certain joints, a command in the form 

ZERO m^, m^, - - 

is followed by a sequence of records in 'loop limit' format (see JREF) 
identifying the joints at which the indicated constraint is imposed. For 
example, to set the direction 3 displacement and direction 2 rotation equal 

to zero at joints 1, 2, 3, 17, 34, 36, and 38, the following language could 

be used. 
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ZERO 3,5: 1,3: 17: 34,38,2$ 


TAB 

CON 


To declare joint motion components which will have non-zero 
values assigned in processor Q, a command in the form 
NONZERO m^, m^, - - - ; 

is followed by input cards in 'loop limit' format identifying the affected 
joints. 


To release components erroneously constrained, the command 
RELEASE m r m^, ~ ; 

is followed by cards identifying the affected joints, as above. 


Additional commands, as summarized below, can be used to 
cause CON to automatically impose constraint to joints lying in global 
PlaneS ' joint in a sym metry or antisymmetry plane must have an 

axis of the joint refer en ce frame (see JREF) perpendicular to the plane , 
although the other two axes need not be parallel to global axes. 


Command 

SYMMETRY PLANE = n 

f 

ANTISYMMETRY PLANE= n 
FIXED PLANE = n 

XZERO= x 


Meaning 


The plane normal to global axis n 
is a symmetry plane. 

The plane normal to global axis n 
is an antisymmetry plane. 

All joints in the plane normal to 
global axis n are completely fixed. 

In processing the above commands, 
x is the linear tolerance used to 
determine if a joint is in a plane. 
Default x is .01. 
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Command 


Meaning 


RZERO= r 


In processing the above commands, 
joints found to be in a symmetry or 
antisymmetry plane must have an 
axis normal to the plane.. The axis 
alignment tolerance is r (default 
value = .0001 radians). 
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3.1.16 JOINT ELIMINATION SEQUENCE (JSEQ) 

As indicated in Reference 1, the sequence in which joints are 
eliminated during the factoring (reduction) sequence performed in 
processor INV significantly affects computer execution time and data 
storage requirements, both in factoring and subsequent execution of 
processor such as DSOL, and EIG. We expect in the future to provide 
software which will substantially automate determination of joint 
elimination sequences. For the present, however, it is necessary for 
the user to do one of the following. 

(1) Carefully study Reference 1 to learn the factors governing 
choice of favorable joint numbering sequences. This is 
not particularly difficult. The rules are less restrictive 
and in some respects simpler than those associated with 
band matrix or wavefront procedures. 

(2) Ask the advice of someone familiar with the sparse matrix 
solution procedure used in SPAR or in the SNAP program. 

For the infrequent user this is an excellent solution, since 
most experienced users can readily identify suitable 
sequences. 

(3) Use a numbering sequence which would be appropriate for 

a band-matrix or wavefront procedure. This will forego the 
opportunity to use the full capability of SPAR's sparse matrix 
procedure, but it is the only alternative to the above. This will 
often give satisfactory costs for problems of small to moderate 
size. 
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JSEQ 


Several nisthods of defining the «ninf a i * * .. 

lining e 3ouit elimination sequence are 

allowed. The simplest and least aittr.rv,=»*^ 

P ana least automatic is comprised of a sequence 

of input cards to string integer format, as illustrated below. 


10/13, 8/5, 9$ 

4, 1/3, 14$ 


or by 


The above implies the following joint sequence: 

10, 11, 12, 13, 8, 7, 6, 5, 9, 4, i >2 , 3 14 

For regular structures involving repeating patterns to the joint sequence, 
the command 

REPEAT n, inc, jmod 

y p cede a packet of cards (terminated by another REPEAT card, 
end of data) defining a list of joints. The effect of the REPEAT card is to 
C ause jmod to be added to each joint to the list. The list is repeated n 
times, With toe added to each joint number upon each repetition. 

For example , suppose it is necessary to describe the following 
lomt elimination sequence for a structure containing 150 joints: 

HI, 112, 113, 120, 

101. 102, 103, 110> 

131, 132, 133, 140, 

121, 122, 123, 23q 

"• 97 — 3, 2, 1. 

The above is represented by the following JSEQ input. 
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. JSEQ 

REPEAT 2 , 20 , 100: 11/20, 1/10: 

REPEAT 0, 0, 0: 99/1$ 

The default values of n, inc, and jmod are 1, 0, and 0, respectively. 
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RMASS 

1* 17 RIGID MASSES (RMASS) 

Two formats are available to define a rigid mass, as indicated below. 

If FORMAT=l, the data sequence on each card is 
k, M, 1^, I^, I^, where 

k is the joint to which the mass is attached, 

M = mass, and 

L = mass moment of inertia about axis i of the joint reference 
frame * jgote that it is requir ed that the joint reference 
frame be coincident with the principal axes of the rigid 
Rigid links and zero-length 2-node elements 
are often useful in modeling large rigid masses. 

If FORMA T= 2, the data sequence is 

I ’ 2 1 ^3/ ^1 1 ^ 2 * ^3* ^here the data has the same where the 

data has the same meaning as above, except that Mi is the mass effective in 
the direction of axis i of the joint reference frame. 


If a data card in either of the above formats is preceded by the 

command 

REPEAT n, j. 

the effect of the REPEAT command is to cause the same mass to be defined 

at n joints, k, (K + j ), <K + 2j), - A REPEAT command applies only to the 

data card immediately following it. 


To cause all subsequently defined masses <M, or Mj, M.,, and M 3 > 
to be multiplied by p, and to cause all I's to be multiplied by q, a command 3 
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RMASS 


is given in the following form. 

CM p, q . 

If successive CM commands are given, the new values of p and q 
replace the previous values. 

The output data set generated by RMASS is, in effect, a diagonal 
mass matrix corresponding to the mass data read from the data cards 
defined above. Distributed mass associated with elements is not 
included in the RMASS output. However, if element distributed mass is 
negligible, RMASS output can be used as a system mass matrix by 
processors Q and EIG. 
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3.2 ELD 

Function - ELD translates element definition data from input cards 
into data sets which are readily usable by other SPAR processors. Elements 
may be defined singly or through a variety of network generators, or com- 
binations thereof. An element is defined by specifying (1) the joints to which 
it is connected, and (2) pointers to applicable entries in tables of section 
properties, material constants, etc. The tables of section properties, etc. , 
are usually generated by the TAB processor. As the ELD input is processed, 
checks are performed to detect errors such as references to nonexistent 
table entries, joint numbers, etc.; however, ELD does not extract any 
data from TAB generated tables . Accordingly, if TAB is subsequently 
re-executed to alter values of section properties, joint position coordinates, 
etc. , it is not necessary to re-execute ELD, provided the connected joint 
numbers, table entry pointers, etc. , of the elements are unchanged. 

Alphanumeric names used to identify specific types of elements 
are summarized in Table ELD-1. Formulations on which these elements 
are based are discussed in Appendix B. 

Table ELD-1: Summary of SPAR Elements 
Description 

General beam elements, such as angles, wide 
flanges, tees, zees, tubes, etc. 

Beams of finite length for which the 6 x 6 
intrinsic stiffness matrix is directly specified. 

Bar elements having only axial stiffness. 


Name 

E21 

E22 
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Name 

Description 

E24 

Plane beam. 

E25 

Zero -length element used to elastically connect 
two coincident joints. 

E31 

Triangular membrane, flat, aeolotropic, using 
TM constant- stress formulation. 

E32 

Triangular bending element, flat aeolotropic. 

E33 

Triangular membrane + bending element, flat, 
aeolotropic, using TM and TPB7 formulations. 

E41 

Quadrilateral membrane, flat, aeolotropic, 
using QMB5 hybrid formulation. 

E42 

Quadrilateral bending element, flat, aeolotropic, 
using QPB11 hybrid formulation. 

E43 

Quadrilateral membrane + bending element, flat, 
aeolotropic, using QMB5 and QPB11 formulations. 

E44 

Quadrilateral shear panel, flat, using QMB1 hybrid 
formulation. 

For purposes 

of explaining card input, we will consider ELD to be 

comprised of an array of sub-processors, one for each type of element. 
The function of each sub -processor is to read input cards defining all of 
the elements of a particular type. A sub-processor is activated by an 
ingut card containing the name of the element type. Sub-nrocessnrs mav 

be called in any orde 

r. A typical input stream is shown below. 

® XQT ELD 


E43 



(Data cards defining all of the E43 elements in the structure) 


E21 

(Data cards defining all of the E21 elements in the structure) 
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Execution of each sub-processor results in production of data sets having 
the following names (EXY represents the name of an element type, such 
as E21, E33, etc. ): 

DEF EXY = Basic (integer form) element definitions. 

GD EXY = Group directory. 

GTIT EXY = Group titles. 

When a sub-processor is executed, the resultant output data sets replace 
all data previously generated by the same sub-processor in any previous 
execution. 

Elements of each type may be assigned by the analyst to separate 
groups . Within each group, each element has an identifying index number . 
The grouping of elements serves various purposes, which new users will 
learn with experience. It is almost always best to use many groups. 
Definition of the elements in group n is preceded by an input card in the 
following form: 

GROUP n' - - - title describing group n - - All elements 
defined by input cards following such a card belong to group n. Groups 
should appear sequentially (group 1, group 2, - - ). If no GROUP card 
appears, all elements are assigned by default to group 1. The order in 
which individual element definitions appear following a GROUP card 
determines the index numbers t hey are assigned within the group. 

* It is planned in a future release to provide an "update" mode of operation 
allowing deletions, additions, and corrections to be made. 
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Element definition cards have the forms shown below. 

All 2 -node elements: 

Jl, J2 

All 3 -node elements : 

Jl , J2, J 3 

All 4-node elements: 

Jl, J2, J3 , J4 

It should be noted that the order in which Jl, JZ - - appear on the element 
. definiti ° n cards is very important, since the orientation of element 
reference frames is defined on the ba sis of this information (see Fig. ELD-1). 
All element-related quantities, such as stiffness coefficients, stresses, etc., 
are relative to these reference frames. 


In the above, Jl, J2, - - indicate the joints interconnected by the first 
element defined by the card. If the optional network operation parameter, 
NETOPT, NET(l), NET(2), - -, are given, they result in definition of 

additional elements, in accordance with specific rules which will be defined 
subsequently. 

The section properties, material properties, etc. associated with 
each element depend upon the values of various table reference pointers 

prevailing at the time the element definition card appears. These pointers 
ere summarized below. 


Optional 


NETOPT, NET(l), NET{2) 


NETOPT, NET(l), NET{2), - - 


NETOPT, NET(l), NET(2), - 
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Table Default 

Pointer Value 

NMAT 1 

NSECT 1 


NOFF 0 

NNSW 0 

NREF 1 


Associated Table 

Material constants, MATC. 

Section properties. E44 elements 
refer to Table SB. All other 3 and 4- 
node elements refer to Table SA. 

E21, E22, E23, E24, and E 25 refer 
to BA, BB, BC, BD, and BB, respectively. 

Beam rigid link offsets, BRL 
(2 -node elements, only). 

Nonstructural distributed weight, NSW. 

Beam reference frame orientation, 

MREF (2-node elements other than E25, 
only). For E25 elements, NREF points 
to an entry in the ALTREF table. 

For 3 and 4-node elements NREF is not 
a table entry pointer. Instead, NREF is 
used to specify the direction of action of 
positive pressure exerted on the element. 
If NREF-O, pressure exerts no force on the 
element. IfNREF=l, positive pressure 
acts in the direction of the 3 -axis of the 
element reference frame (see Fig. ELD-1, 
below). If NREF=-1, pressure acts in the 
opposite direction. 



(1) All systems right-hand 
rectangular. 

(2) Orientation of 1, 2 axes of 
2 node elements depends on the 
entry in MREF table indicated by 
NREF. 

For 3, 4-node elements, J3 lies 
in quadrant 1 of plane 1-2. 


3 



Figure ELD-1. Element Reference Frames 
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An example of the use of the table reference parameters is shown 
below. 



NMA T= 4: NSECT-7: 29,13: 13,74: NSECT=8: 74, 16: NMAT=2: 74,92$ 


ELD also accepts input commands which have the effect of modifying 

the data given on subsequent input cards. These commands are summarized 
below. 


Command 


Meaning of n ( n default is zero in all cases) 


MOD 

MOD 


MOD 

MOD 

MOD 

MOD 

MOD 

INC 

INC 

INC 

INC 

INC 


JOINT = 

n 

Add 

n 

to all joint numbers. 

GRO UP= 

n 

Add 

n 

to all group numbers, 

NSECT= 

n 

Add 

n 

to NSECT pointers. 

NMAT= 

n 


1 1 

NMAT " 

NNS W= 

n 


1 1 

NNSW 

NREF= 

n 


n 

NREF » 

NOFF= 

n 


1 1 

NOFF " 

NSECT= 

n 

See below. 

NMAT= 

n 

11 

1 1 


NNSW= 

n 

1 I 

if 


NREF= 

n 

ri 

it 

* 

NOFF= 

n 

1 1 

1 1 

* 
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A frequently encountered situation is the one in which a substantial 
portion of a finite element model is comprised of a regular mesh, which 
can be automatically defined by a mesh generation command, but the 
section properties are so nonuniform that it becomes necessary or highly 
desirable for simplicity of input preparation to associate with each element 
a unique entry in the section property table. A common example is an air- 
craft fuselage, where cutouts, doublers, stiffeners, etc., result in extensive 
variations in section properties over a mesh that is topologically and geo- 
metrically regular. A similar situation involving distributed nonstructural 
weight also sometimes occurs; e.g. thermal protection material of varying 
thickness. The INC command is intended primarily to address situations 
of this type. 

The INC commands are used to cause the associated table reference 
pointers to be automatically incremented by n as successive elements are 
defined. This is especially useful in conjunction with network generation 
options, which are defined in detail on subsequent pages. For example, 
the element definition command 701, 702, 1 , 50 causes fifty elements to be 

defined, connecting joint pairs (701, 702), (702,703) (750,751). If 

section property table entries 101, 102, 150 apply, successively, to 

these elements, the following input would s uff ice. 

NSECT= 101: INC NSECT=1: 701,702, 1,50$ 

All MOD and INC parameters are internally reset to zero after 
conclusion of each sub-processor execution . The MOD and INC commands 

are not cumulative. That is, MOD NSECT 7: MOD NSECT 3: is 

not equivalent to MOD NSECT 10. 
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Reset Controls 

Only one reset command is permitted: 

OUTLIB = Destination library for ELD output. The TAB-generated 
data sets must already be resident in this library. 

Core Requirements 

ELD rarely requires more than 1500-3000 words of working core 
space, regardless of the size of the problem. 

Code Release Data 

Level 7 (UlNttVAC, CDC), July 1974, coded by W. D. Whetstone. 
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Two-node element network generators . 


If NETOPT=l, . > 

NI= NET{I) (defaults), 

NJ= NET{2) (default=l), 

JINC= NET<3). 

Implied sequence: 

N1=J1 

IDIFF= J2-J1 
DO 200 J= 1 , NJ 
DO 100 1=1, NI 

N2=N1+IDIFF 

Define element connecting node Nl to N2 
N1=N2 

N1=J1 + J * JINC 
Example: 


100 

200 


J1 

9 


32 

12 


Netopt NI NJ JINC 

1 4 2 100 


9 (1)* 12 


( 2 ) 




21 



ELD 



109 112 

' (67 


115 118 121 

~*~Tn ' ( 8 ) 


‘'The order in which the elements are defined is indicated h v „ k 
enclosed in parentheses. The index number identifying element^^hhi 
each group are determined by the order in which the^fements are defSe, 
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Two-node element network generators (continued). 


If NETOPT=2 , 

NI= NET(l) (defaults 1), 

NJ= NET(2) (defaults 1 ) , 

J1NC= NET(3). 

Implied sequence: 

Nl= J1 
IDIFF= J2-J1 
DO 200 J =1 , NJ 
DO 100 1=1, NI 

N2 = NI + IDIFF (except for closing element, when I=N1) 
Define element connecting N1,N2. 

100 N1=N2 

200 N1=J1 + J*JINC 


Example: 






J1 

J2 

Netopt 

NI 

NJ 

JINC 

9 

12 

2 

4 

2 

100 
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Two-node element network generators (continued). 


If NETOPT=3 , 

NI= NET(l), 

UNC= NET<2), 

NJ= NET(3) < defaults) 
JINC= NET{4). 


Implied sequence: 


N1=J1 
N2= J2 

DO 200 J=l, NJ 

DO 100 1=1, NI 

Define element connecting Nl to N2 
100 N2 = N2 + IINC 

Nl = J1 + JINC 

200 N2 = J2 + JINC 


Example: 

J1 

J2 

Netopt 

NI 

IINC 

NJ 

JINC 

4 

37 

3 

4 

10 

2 

100 
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Three-node element network generators. 


If Netopt= 1 , 

NI= NET(l), 

NJ= NET(2) {default= 1 ). 

Implied sequence: 

HNC= J2-J1 
JINC= J3-J2 
NI= J1 
N2= J2 
N3= J 3 
N4= J3-IINC 
DO 200 J=l, NJ 
DO 100 1=1, NI 

Define element connecting Nl, N2, N3 
Define element connecting N3, N4, Nl 
Nl = Nl + HNC 
N2= N2 + IINC 
N3= N3 + XINC 
100 N4= N4 + IINC 

Nl= J1 + JINC 
N2= J2 + JINC 
N3= J3 + JINC 
200 N4= N3 - IINC 


Example: 


J1 J2 J 3 


Netopt NI NJ 


2 3 103 


3 2 
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Three-node element network generators (continued}. 



NI= NET(l), 

NJ= NET(2) (default=l), 
JINC= NET(3). 


Implied sequence: 


HNC= J3-J2 
Nl = J1 
N2= J2 

DO 200 J=l, NJ 

DO 100 1=1, NX 

N3=. N2 + IINC (except closure when I=NX) 

Define element connecting N1.N2.N3 

100 N2=N3 

N1 = J1 + J*JINC 
200 N2 = J2 + J*JINC 


Cv 


Example: 

J1 J2 J 3 
2 5 7 


Netopt NX NJ 

2 6 Z 


JINC 

30 
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Four -node element network generators ♦ 


If Netopt= 1 , 

NI= NET{1), 

NJ= NET(2) (default^!), 

NK= NET{3) (default I), 

KINC= NET (4). 

Implied sequence : 


IINC= J2-J1 
JINC= J4-J1 
DO 300 K=l, NK 
N1 = J1 

DO 200 J = l, NJ 
DO 100 1=1, NI 

N2= Nl+IINC 
N3= N2+JINC 
N4= Nl+JINC 

Define element connecting N1,N2,N3,N4 
100 Nl = Nl+IINC 
200 Nl = J1+J*JINC 
300 J1 = Jl+KINC 

Example* 1 ': 


J1 

J2 

J3 

J4 

2 

3 

23 

22 

62 

63 


64 


42 

(5) 

43 

(6) 

44 


(3) 


(4) 

22 


23 

24 j 

(1) 

(2) 


3 4 


NETOPT NI NJ 

1 2 3 

262 263 


NK KINC 
2 200 

264 


242 

222 


202 


(11 > 243 

(12) 

223 j 

(10) 

i 

(7) | 

(8) 


244 

224 


203 


204 


Note: J3 must be present, although not used. 
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Four -node element network generators (continued). 


If Netopt=2 , 

NI= NET(l), 

NJ= NET(2) (default=l). 

Implied sequence: 

XINC= J2-J1 
JINC= J4-J1 
Nl= J1 
DO 200 J = 1 , NJ 
DO 100 1=1, NI 

N2 = N1 + UNO 

N3 = N2 + JINC (except closure when I=NX) 
N4 = Nl + JINC (except closure when I=NI) 
Define element connecting N1,N2,N3,N4 
100 Nl = Nl + IINC 
200 Nl = J1 + J*JINC 


Example: 


J1 J2 J3 J 4 Netopt NI NJ 
1 11 12 2 2 63 


3 13 23 


33 



32 

31 
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3.3 TO£ O. 

Function. TOPO analyzes element interconnection topology and creates 
two data sets, KMAP and AMAP, KMAP is used by processors K, M 
and KG to guide assembly of system stiffness and mass matrices in the 
SPAR standard sparse -matrix format. AMAP is used by processor INV 
in factoring system matrices, input data sets are: 

# All basic element definition data sets generated by 
processor ELD. 

• JSEQ, the joint elimination sequence, generated by TAB 
sub-processor JSEQ. 


RESET Controls: 


Name 

D efault 
Value 

Meaning 

BLIB 

1 

Library containing all input data sets 

HLIB 

1 

Destination library for KMAP 

ILIB 

1 

Destination library for AMAP 

LRKMAP 

1344 

KMAP block length 

LRAMAP 

1344 

AMAP block length 

FRTKMAP 

0 

Nonzero value results in printout of 
topological data during generation of 
KMAP. 

PRTAMAP 

0 

Same as above for AMAP. 

MAXSUB 

1000 

Maximum allowable value of size index 
I g (see Ref. 1). 

ILMAX 

(none) 

Maximum joint connectivity to be permitted 
in allocating space for topological tables. 
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Name 


Default 

Value 



TIME 0 


Nonzero value causes printout of intermediate 
CP and wall clock times . 


Core re quirements . Core space requirements for TOPO are not 
easily defined because they depend on the number of nodes, element 
interconnection topology, and the joint numbering sequence prescribed 
The working core required is not often greater than 

3000 + 2 * J + 4 * 1^, where 

J = the number of joints, and 
I = maximum joint connectivity. 

Code release data. Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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Function. For each element in the structure there is an "element 
information packet. " Depending upon the particular type of element, 
each packet usually contains information in each of the following 
categories: 

(1) General information such as the connected joint numbers, 
table reference integers, etc. 

(2) Material constants. 

(3) Geometrical data, e.g. dimensions, transformation 
matrices, etc. 

{4) Section property data. 

(5) Intrinsic stiffness matrix (or equivalent). 

(6) Intrinsic stress matrix (or equivalent). 

Processor E constructs these packets in skeletal form (omitting the 
intrinsic stiffness, and stress data). Input to £ consists of* 

• All TAB-produced data sets except for rigid mass data, 
constraint data, and the joint numbering sequence. 

• All data sets produced by ELD. 

Output from E consists of two data sets for each element type present 

Where Exx represents the element type (e.g., E21, E22, E31, 

E41, E42, etc.), the output data sets are 



E 

(1) Exx = all element information packets for 

element type Exx, and 

(2) DIR Exx = an associated directory. 
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RESET Controls. 


Name 

Default 

Value 

Meaning 

BLIB 

1 

Library containing all input data sets. 

ELIB 

1 

Destination library for all output data sets. 

LZERO 

. 001 

Used to detect identically zero lengths. 

WARP 

.001 

Test for allowable warp of quadrilateral 
elements, such as E41, which are required 
to be flat. Warp is defined as Z/L, where Z 
is the distance from the fourth node to the 
plane of the first three nodes, and L is the 
distance between nodes 1 and 2 . 

RCH 

. 0001 

Where R is any 3x3 orthogonal trans- 
formation, terms of R R must equal the 
identity matrix, within a tolerance of RCH. 

TIME 

0 

Nonzero value causes printout of inter- 
mediate CP and wall clock times. 

Core Requirements, 

It is rarely necessary to allocate more than about 


20, 000 words of working core space for execution of E. Requirements 
vary from problem to problem depending upon the size of the TAB- 
generated tables. All of the geometry, material constant and distributed 
weight tables are held continuously in core during execution of E. As 
packets are formed for each successive element type, the associated 
section property table is loaded into core. In addition to the above, 
about 1000 words should be allowed for input of the basic element 
definition data sets from ELD, plus 5000-10000 words to serve as an 
output buffer for the element information packets. 

C ode release data . Level 7 (Uni vac, CDC), July 1974, coded by 
W. D. Whetstone. 
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3.5 EKS. 

Function. Based on the dimensions, section properties, etc., currently 
imbedded in the element information packets originated by processor E, 
EKS computes intrinsic stiffness and stress matrices and inserts them 
into the packets. 

RESET Controls. 

Default 


Name 

Value 

Meaning 

ELIB 

1 

Library containing the element information 
packets. 

TIME 

0 

Nonzero value causes printout of intermediate 
CP and wall clock times. 

gazero 

10. - 20 

Zero-test parameter, (beam area) x (shear 
modulus ) . 

Cl ZERO 

10. - 20 

Zero-test parameter , beam non-uniform 
torsion constant. 


Core Require ments . EKS requires only a buffer area through which 
element information packets are transmitted. About 5000-15000 
locations are usually suitable. IO counts will vary in inverse proportion 
to core space. 

Code release data . Level 7 (Univac, CDC), July 1974, routines 

QMB1, TM, TPB7, QMB5, QPB11 coded by C. L. Yen, all others 
coded by W. D. Whetstone. 


3. 5-1 



K 


3. 6 K_ 

Function. K assembles unconstrained system stiffness matrices 
in the standard SPAR sparse-matrix format. Input consists of the 
following: 

9 The element information packets produced by processors E 
and EKS. 

® Topological information contained in data set KMAP 
produced by processor TOPO. 

• Certain minor directory information generated by 
processors TAB and ELD. 

The output is a data set named K SPAR, the unconstrained system 
stiffness matrix. 


RESET Controls. 


Name 

BLIB 

ELIB 

HLIB 

KLIB 

LREC 


Default 

Value Meaning 

1 Library containing TAB, ELD output. 

1 Library containing element information 

packets. 

1 Library containing data set KMAP, produced 

by TOPO. 

1 Destination library for K output data. 

2240 Block length of output data set. 


Core Requirements . It is rarely necessary to allocate more than 
10000-20000 locations. Space is req\iired to concurrently accommodate 
the following: 
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• One block of AMAP 

• One block of K SPAR 

• A core pool area in which intermediate results are held 
during assembly of the system matrix. Specific require- 
ments depend on element interconnection topology and the 
joint numbering sequence. Required space will be much 
less than I g submatrix storage blocks (see Ref. 1). 

Code Release Data . Level 7 (Univac, CDC), July 1974, coded 
by W. D. Whetstone. 
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3. 7 M_ 

Function. M assembles unconstrained system consistent mass 
matrices, considering only the structural and nonstructural distributed 
mass associated with the elements. That is, rigid masses lumped at 
joints are not included in M SPAR, the output data set produced by 
processor M. The first production release of M (level 7) ignores rigid 
link offsets in its calculations. 

M requires the same input data sets as K, and has identical 
working core requirements. 


RESET Controls. Reset controls are the same as for processor 
K, plus the following: 


Name 

G 


IBEAM 


INERT 


Default 

Value Meaning 

1. Constant used to convert element distributed 

"weight" to mass. For example; if distributed 
weights were given in lbs/in and lbs /in 2 , and 
material specific weights were given in lbs/in^. 
one would normally reset G = 386. , depending 
upon the planned use of M SPAR. 

0 Any nonzero value will cause all beam rotatory 

inertia terms, both lateral and torsional, to 
be ignored. 

0 Any nonzero value results in assumption of 

linearly varying lateral displacement fields 
in calculating mass matrices for all elements. 
That is, beams are treated as axial (bar) 
elements, shell elements are treated as 
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Default 

Name Value Meaning 

membranes, etc. This often is an excellent 
approximation, and requires only a small 
fraction of the CP time used in the default 
mode. 


Code Release Data. Level 7 (Uni vac, CDC), July 1974, routines 
M33, M34, M63 and M64 coded by C. L. Yen, all others coded by 
W. D. Whetstone. 
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3. 8 KG 

Function. KG assembles unconstrained system initial- stress 
(geometric) stiffness matrices, K , in the standard SPAR sparse matrix 

O 

format. The first production release of KG (level 7) ignores rigid link 
offsets in its calculations. The output data set is named KG SPAR. Required 
input data sets are the same as for processor K, plus one additional data set 
containing the pre -stress solution upon which K is to be based. Core 
requirements are the same as for K, plus additional space (not more than 
7 times the number of joints) in which the pre- stress solution is resident. 

RESET Controls. Reset controls are the same as for processor K, 
plus the following: 

Default 


Name 

Value 

Me aning 

QLIB 

1 

Library containing the pre- stress solution 
on which Kg is to be based. 

QSEQ 

1 

Sequence number identifying the solution 
case on which Kg is to be based. 

IKG2 

0 

Any nonzero value results in all beams 
being treated as though they were axial 
elements (bars) in computing Kg. 

IKG34 

0 

-10 

Any nonzero value results in all shell 
elements being treated as though they were 
membranes in computing Kg. 

ZERO I 

10. 

Zero-test parameter used to detect identically- 
zero beam moments of inertia. 

zeroa 

10. " 10 

Zero-test parameter used to detect identically - 
zero cross-sectional areas of beams. 

Code release data. 

Level 7 (Univac, CDC), July 1974, routines 


GTM, GQM, GTP, GQP coded by C. L. Yen, all others coded by 
W. D. Whetstone. 
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3. 9 LCM 

Function.. LCM forms linear combinations of system matrices 
in the SPAR standard sparse matrix format. That is, where Q1 and Q2 

are existing system matrices, and Cl and C2 and constants, LCM forms 
the data set R: 

R = C1*Q1 + C2*Q2. 

The default operation of LCM is to add K to Kg. 


RESET Controls . 

Default 


Name 

Value 

Q1LIB 

1 

Q2LIB 

1 

RLIB 

1 

Cl 

1.0 

C2 

1. 0 

Ql 

K 

Q2 

KG 

R 

K+KG 


Meaning 

Library containing Ql. 
Library containing Q2. 
Destination library for R. 
Value of Cl 
Value of C2 

First word in name of Ql. 
First word in name of Q2. 
First word in name of R. 


Example. Calculate 3. 7 M + 20. OK, where M SPAR is in 
library 5 and K SPAR is in library 7. The output is to be named CMK 
SPAR, and stored in library 2, 

@XQT LCM 

RESET Cl = 3. 7, Ql = M, C2 = 20. , Q2 = K 
RESET R = CMK 

RESET Q1LIB = 5, Q2LIB = 7, RLIB = 2 
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Core Requirements . The working core needed is 3 times the 
block length of Ql. Q2's block length must be equal to that of Ql. 

Code Release Data. Level 7 {Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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3. 10 MRM 

Function. The function of MRM is to add a diagonal matrix to a 
system matrix in SPAR's standard sparse matrix format, resulting in the 
production of two data sets (1) the sum of the two matrices in SPAR format, 
and (2) the diagonal of the sum of the two matrices. 

The normal application of MRM is to add a diagonal mass matrix 
named RMASS (produced by TAB sub-processor RMASS) to the assembled 
system consistent mass matrix M SPAR, naming the sum M+RM SPAR. 
The data set containing the diagonal of the sum is M+RM DIAG. 


RESET Controls. 


Name 

Default 

Value 

Meaning 

MLIB 

1 

Library containing M SPAR, 

RMLIB 

1 

Library containing RMASS. 

LIB 

1 

Destination library, output data sets. 


Core Requirements . Where J is the number of joints, and L is 
the block length of M SPAR, the core required is J + L. 

Code Release Data. Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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3.11 inv 

Function. INV factors assembled system matrices in the SPAR 
standard sparse matrix format, using a procedure essentially the same 
in principal as the one outlined in Ref. 1. The following input data sets 
are required. 

• An unconstrained assembled system matrix, such as 
K SPAR. 

• A constraint definition, such as CON 1, produced by the 
CON sub-processor of TAB. 

• Topological information contained in data set AMAP 
produced by processor TOPO. 

Output consists of a factored system matrix having the name INV XXXX n, 
where 

XXXX = the first word of the name of the input stiffness matrix, and 
n " an integer identifying the constraint case. 

For example, if K SPAR is factored using constraint case 1, the 
output is INV K 1 . 


RESET Controls 


Name 

Default 

Value 

Meaning 

KLIB 

1 

Library containing input system matrix. 

KILIB 

1 

Destination library for INV output. 

ILIB 

1 

Library containing topological data set AMAP. 
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Default 
Name Value 

K K 

CON 1 

LRA 3136 

DZERO .00001 


Meaning 


First word in name of input data set. 
Typical RESET actions are K = K+KG 
K = K+KA, K = K+CM, etc. 


Integer identifying the constraint case 
to be used. This data set must be 
resident in library 1. 

Output data set block length. 

Zero-test parameter used to detect 
identically zero diagonal terms during 
factoring of the system matrix. 


Core Requirements . 

Where J = the number of joints in the structure, 

L.^ = block length of input matrix, 

= block length of output matrix, 

= block length of AMAP and 

no. of degrees of freedom per joint (usually 3, 4, 
5 or 6), and 


n = 


- the size index (Ref. 1) of the structure, the minimum 
core required is 

J + + L.^ + L^ + n^I g . 


Code Release Data. Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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Function. Q reads card input defining static load cases. The 
card data is translated into various data sets which are subsequently 
stored in library 1, or any other library selected by the user via RESET 
Control. The corresponding static solutions are subsequently computed 
by other processors. The definition of a static loading case is initiated 
by a command in the following form: 

CASE n' ALPHA -NUMERIC TITLE f 

where n is an integer uniquely identifying the load CASE. It is advisable, 
though not required, for CASE 1 to appear first, CASE 2 second, etc. 
Following the CASE card there are one or more packets of cards, each 
packet defining one of several types of loading (e. g. pressure, thermal, 
etc. ) which together comprise the CASE. A packet is initiated by any of 
the following multi-word commands: 

NODAL FORCES 

NODAL MOTIONS 

NODAL TEMPERATURES 

NODAL PRESSURES 

A particular type of packet may not appear more than once within 
each CASE. There is no limit on the number of data cards in each packet. 
The order of individual data items on each card is defined below. Omitted 
trailing items are assumed zero. 

NODAL FORCES: 

jo, k, q, ni, iinc, dqi, nj, jinc, dqj. 

Interpretation: If k = 1, 2, or 3, q is a concentrated force 

acting on joint jo in direction k, relative to the joint reference 
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frame associated with joint jo. If k = 4, 5, or 6, q is a 
moment in direction k - 3. 

If given, ni, iinc, etc. , imply the following: 

For J = 1 through nj and 1=1 through ni, a direction 
k load of magnitude F is applied to joint N, where 
N = jo + (J - 1) * jinc + (I - I) * iinc, and 
F = q + (J - 1)* dqj + (!-!)❖ dqi. 

Note that the orientation of individual joint reference frames 
is prescribed via the TAB processor. 

NODAL MOTIONS: 

Card input is exactly the same as for NODAL FORCES, except 
that q is interpreted as a specified nodal displacement or 
rotation. Only motion components declared to be NONZERO in 
the CONSTRAINT subprocessor of TAB may have non-zero 
prescribed values. Any motion components declared NONZERO 
in TAB, but not assigned a value for a particular CASE, are 
set equal to zero. 

NODAL TEMPERATURES: 

jo, q, ni, iinc, dqi, nj, jinc, dqj. 

Interpretation: q is the temperature of joint jo. Parameters ni, 
iinc, etc. serve the same function as in the NODAL FORCE input. 
Note that through RESET control, temperatures may be shifted 
and scaled. 
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NODAL PRESSURES: 

Data card interpretation is the same as for NODAL TEMPERATURES, 
except that q is nodal pressure. The positive direction of pressure 
acting on a 3 or 4 node element is determined by the value of NREF 
assigned when the element was defined in ELD. If NREF is zero, 
no pressure acts on the element. If NREF = 1, positive pressure 
acts in direction 3 of the local element reference frame. If 
NREF = -1, positive pressure acts in the opposite direction. 

Pressure scaling and shifting are available through RESET controls. 

In addition to the above, each CASE may include inertia loading 
defined by one or more commands of the following form: 

INERTIA FORCE, k, acc. 

The above indicates linear acceleration, acc (not a "g" factor), 
in global direction k (k = 1, 2, or 3). Note that k must be an 
integer, and acc must be real. 


Inertia loads are formed by linearly combining inertia load 
influence vectors named RBMX VEC (n = 1, 2, 3). If these 
vectors do not already reside in the Q output library, it will be 
necessary to RESET MLIB = k, where k is the library containing 
the system mass matrix to be used in calculating the inertia load 
influence vectors. In the default mode, inertia forces are based 
on the assembled system mass matrix M SPAR (see Processor M). 
Through RESET control other choices are available. 
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In future releases, a variety of additional methods are planned 
for defining all of the above, and other types of loading. 

n, Q produces the following data sets. 
Contents 

Alphanumeric title supplied on the input 
card initiating CASE n. 

Nodal temperatures (if given). 

Nodal Motions (if given). 

Total equivalent nodal forces and moments 
from all sources. 

In addition, a data set named LDIR is produced, containing the 
smallest and largest CASE number processed in the execution. LDIR 
is used in the default operation of several processors. 

RESET Controls. 

Default 


Name 

Value 

Meaning 

ELIB 

1 

Library containing the element information 
packets generated by processors E and EKS. 

MLIB 

0 

Library containing the system mass matrix 
to be used in calculating inertia forces. 

QLIB 

1 

Destination library for all output data sets. 

TMULT 

1.0 

Input nodal temperatures are scaled 
(multiplied) by TMULT. 

TNULL 

. 0 

TNULL is subtracted from scaled nodal 
temperatures. 

PMULT 

1.0 

Input nodal pressure are scaled (multiplied) 
by PMULT. 


For load CASE 

Name 
CASE -- n 

NTMP -- n 
NMOT -- n 
NFM - - n 
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Name 

Default 

Value 

PNULL 

. 0 

M 

M 

MTYP 

SPAR 


Core Requirements. 


Meaning 


PNULL is subtracted from s called 
nodal pressures. 

See below. 

The default name of the mass matrix 
input is M SPAR. For example, if 
you elect to use consistent mass 
matrix M+RM, RESET M = M+RM. 

If you elect to use only the rigid mass 
data input in TAB, RESET M = RMASS, 
MTYP = MASK. 

Working core requirements normally will 


not substantially exceed 12 to 15 times the number of joints in the structure. 


Code Rele ase Data . Level 7 (Univac, CDC), July 1974, coded 
by W. D. Whetstone. 
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3.13 DSOL 

Function . DSOL computes static solutions in the form of joint 
motion vectors. Input to DSOL consists of the following. 

• An assembled system matrix (default is K SPAR), and its 
pseudo -inverse (default is INV K 1), calculated by INV. 

• Data sets comprised of equivalent nodal forces and moments, 
and/or specified nodal motions. These are normally data 
sets produced by processor Q, bearing the names NFM --n, 
and NMOT -- n, where n is the load CASE identifier. 

• The data set LDIR produced by Processor Q, identifying the 
smallest and largest CASE numbers processed in the last 
execution of Q. 

The joint motion data output corresponding to CASE n is named SSOL -- n. 

RESET Controls. 

Default 


Name 

Value 

Meaning 

KLIB 

1 

Library containing systems stiffness matrix. 

KILIB 

1 

Library containing factored system stiffness 
matrix. 

QLIB 

1 

Library containing loading data produced by 
processor Q, into which DSOL will insert 
its results. 

K 

K 

First word of the name of the system 
stiffness matrix. 

CON 

1 

Constraint case used in factoring the system 
stiffness matrix. 

LI 

None 

See below. 
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Default 

Name Value 

L2 None 


ROUNDOFF .0001 


MAXIT 1 


EP 1 

UP 0 

FP 0 

TIME 1 


If both LI and L2 are given, load CASES 
LI through L2 are processed. If LI 
alone is given, only that case is processed. 
If neither LI nor L2 is given, all CASES 
generated in the last execution of Q are 
processed. 

Maximum normalized error to be 
permitted in displacement solutions. 

Data sets corresponding to non-qualifying 
solutions are markedin error. 

Maximum number of iterations allowed 
in attempting to recover from round-off 
error. A value of zero causes nodal 
reaction - force check calculations to 
be omitted. 

Nonzero value results in printout of 
normalized errors. 

Nonzero value results in printout of 
joint motions. 

Nonzero value results in printout of joint 
force summations (reactions, nodal force 
balance errors). 

Nonzero value causes printout of 
intermediate CP and wall clock times. 


Core Requirements. Working core required will not significantly 

exceed 12 times the number of joints, plus the block length of the factored 
system stiffness matrix. 


Code Release Data. Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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3. 14 LCD 

Function. The primary function of LCD is to form linear 
combinations of solutions in the form produced by processors Q and 
DSOL. Card input to LCD consists of a sequence of packets. Each 
packet may be in the following form. 

CASE j' -- alphanumeric title for case j 

N r Cj 

L 2- N 2- C 2 


The above is interpreted as follows: Form a new solution case, identified 
by the integer j, comprised of 

C. times solution case N from library L 

1 1 

+C 2 times solution case N ? from library L 
+ 

+C r times solution case from library L . 

If C is omitted, it is assumed to be 1. 0. The alphanumeric label given 
on the first card of the packet serves the same purpose as in processor Q; 

i. e. to identify CASE j. LCD can be used to combine results produced by 
LCD. 


The following more general packet form is also permitted. 
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CASES 


L 1 N 1 
L Z N 2 


L N 
r r 


j, k 1 title {applies to cases 



C 

r 


I 

I 


1 

2 


I 

r 


j 


through k ) 


The above language results in generation of a succession of combined 

cases: j, j + 1, k. Case j is constructed as previously defined. 

Then I. is added to N. (for i = 1, 2, ---r), and case j + 1 is formed by 

the same rule, etc. The default value of I. is 1. Each case, j, j + 1 k, 

is given the same title from the label field of the CASE card. 

For example, suppose it is desired to define combined cases 40, 41, 
47 as follows: 

Case 40 - 4. 5 (case 130) +2.1 (case 22) + 6.7 (case 204) 

Case 41 = 4. 5 (case 131) + 2.1 (case 24) + 6.7 (case 204) 


Case 47 4. 5 (case 137) + 2.1 (case 36) +6.7 (case 204) 

If cases 130 - 137 are in library 3, cases 22 - 36 in library 2, and case 
204 in library 5, the following packet produces the desired result. 

CASES 40, 47' TITLE 

3 130 4.5 

2 22 2. 1 2 

5 204 6. 7 0 
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The above language is useful in executing extensive sequences of 

combined loads (e.g. symmetric + antisymmetric, thermal + mechanical, 
static + dynamic, etc. ) 


The title of existing case n, originated either in Q and DSOL or 
in LCD, may be changed by a card in the following form. 


CASE n' New Title 


RESET Controls. 


Name 


Default 

Value 


Meaning 


OUTLIB 

THERM 


Output data set destination library. 

If no thermal loading is present, setting 

ITHERM — 0 will slightly reduce execution 
costs. 


Core Requirements . Working core requirements will not significantly, 
exceed 13 times the number of joints in the structure. 

Code ReleaseData. Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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3. 15 GSF 

Function. GSF generates element stresses. Input to GSF normally 
consists of the data sets generated by Q and DSOL, or combined solutions 
produced by LCD. GSF also accesses the element information packets 
originated by processors E and EKS. Output data sets have names of the 
form STRS EXXX N, where EXXX identifies the element type, and N 
identifies the solution case. If no input cards follow the @XQT GSF 
command (and any RESET cards), stresses are computed for all elements 
in the structure. If it is desired to limit stress calculation to designated 
types of elements, one or more input cards in the following form must 
appear. 


ELTS EXXj EXX Z 

where EXX. is the name of an element type for which stress data is to be 
computed. 

RESET Control. 

Default 


Name 

Value 

Meaning 

ELIB 

1 

Library containing element information 
packets, etc (from ELD, E, EKS). 

QLIB 

1 

Library containing solution cases for 
which stress data is to be computed 
and into which stresses are to be stored. 

LI 

None 

Same as in DSOL. 

L2 

None 

Same as in DSOL. 

LREC 

5600 

Block length of output data sets. 
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Name 

Default 

Value 

Meaning 

TIME 

0 

Nonzero value causes printout of 
intermediate CP and wall clock time. 


Core Requirements . Working core requirements will not significantly 
exceed (1) seven times the number of joints plus (2) the block length of the 
output data (default = 5600), plus (3) a buffer area for handling the element 
information packets. 5000-10000 is usually a suitable size for the buffer. 

Code Release Data . Level 7 (Univac, CDC), July 1974, routines 
PBSR3, PMSR4, PBSR4 coded by C. L. Yen, all others coded by 
W. D. Whetstone. 
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3. 16 PSF 

F m ic tion. PSF prints stress data generated by GSF. Input normally 
consists of the data sets generated by Q, GSF, and ELD. If no input cards 
{other than RESET controls, if any) are given PSF will print stress data 
for all elements. If it is desired to limit stress print to designated element 
types, one or more input cards in the following form must appear. 

ELTS EXXj EXX 2 

where EXX. is the name of an element type for which stress data is to be 
printed. 


RESET Controls . 


Name 

Default 

Value 

Meaning 

QLIB 

1 

Library containing output from Q and 
DSOL. 

LI 

None 

Same as DSOL and GSF. 

L2 

None 

Same as DSOL and GSF. 

LINES 

56 

Lines per printed page. 

IE A 

1 

Zero value will cause an abort if input 
errors are detected. 

format 

1 

Print format controller. A value of 2 
results in short-form printout for certain 
of elements. 


Cora Requirements. Working core space does not substantially 
exceed the block length of GSF output. 

Coda Release Data. Level 7 (Univac, CDC), July 1974, coded by 
W . D. Whetstone. 
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3.17 EIG 


F unction. EIG solves linear vibration and bifurcation buckling 
exgenproblems of the types indicated by eqs. (1) and (2). 


rMX - KX = 0 

( 1 ) 

rKgX + KX = 0 

( 2 ) 

K and Kg must be in the SPAR standard sparse matrix format. M may 
be m either SPAR sparse matrix or diagonal format. Although K must 
be nonsingular, it need not be positive definite. 


EIG implements an iterative process consisting of a Stodola 
(matrix iteration) procedure followed by a Rayleigh-Ritz procedure, followed 

by a second Stodola procedure, etc., resulting in successively refined 

approximations of m eigenvectors associated with the m eigenvalues of 

eqs. (Dor (2) closest to aero. Clos ely spaced roots do not adversely 
affect the process* 


In the following discussion of application of the process to eq. (1), 

it is assumed that M and K are n x n, and that m iinearly independent ’ 
system veetors, Y*. Y 2 , . . y™, are Wn _ Meth(jds q( ^ 

vectors will be discussed later. Usually m is chosen to be much less 

than n; that is, m is usualiy 4 to 30, while n may be extremely large 
(e. g. 10000 +). 


In the following discussion, Z is a general linear 
of the Y's: 


combination 


Z = q l Y 1 + <>2 Y 2 + 


+ q Y 

m m 


= HQ, where 
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H = (Y. Y- Y ), and 

1 2 m 

Q = (q t q 2 — - q m )*. (3) 

The Rayleigh -Ritz procedure consists of replacing X with Z in eq. (1); 
that is, substituting (3) into (1) and pre -multiplying by H*, 

r (H*MH)Q + (H*KH)Q = 0 ( 4 ) 

Using the Cholesky-Householder method, this low-order eigenproblem 

is solved for all m eigenvectors Q 1> Q 2> Q^, The Rayleigh -Ritz 

procedure is concluded by using eq. (3) to calculate m improved 
approximations of the n - order system eigenvectors. 



Z = 
m 


HQ 


m 


(5) 


The Stodola (matrix iteration) step is as follows,. From each Z 
of eq.. (5) a new Y is computed, subject to the requirement that 

MZ - KY. (6) 

In performing MZ calculations and in solving for Y, SPAR's sparse- 
matrix algorithms are used. The iterative process continues by using 
the new Y's in another Rayleigh-Ritz procedure, etc; Vectors are 
regularly renormalized to avoid scaling problems. 

The convergence resulting from eq. (6) is readily observed by 

considering vectors Z and Y as linear combinations of the n system 

eigenvectors X., X, -- X . 

I 2 n 
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Z 


n 



a. X. , 
1 i 


n 




( 7 ) 


Substitution of (7) into (6) and 


n 


£ 

i=l 


X* M X. a. 
J 1 1 


pre -multiplication by X? gives 



X* K X. b.. 
J i i 


Since, for i not equal to j, 


X* M X. = X* K X. = 0, 
J 1 3 i 


X* M X. 

b. = -4 i 

X . K X. 

i l 


a. 

i 


r. 

l 


(8) 


where r. is the eigenvalue associated withX^ Eq. (8) indicates that in 
the transition from Z to Y the magnitudes of corresponding eigenvector 
coefficients are diminished in inverse proportion to the associated eigen- 
value magnitudes. Suppose, for example, that some of the system eigen- 
values are r ^ = . 5, r ? = 10. r 25 - 1000. Then 

b, = 2.0 a., 

1 1 

. 1 3Lj, and 
b 25 = *° 01 a 25' 


Since K need not be positive definite, EIG can be used to solve 
"shifted" eigenproblems . For example, if cMX is added and subtracted 
from the left side of Eq. (1), 

(r - c) MX - (K - cM) X = 0. ( 9 ) 

Eq. (9) is of the same form as eq. (1), except that (r - c) has replaced 
r, and (K - cM) has replaced K. Shifting has several useful applications. 
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For example, suppose it is desired to compute only the vibrational modes 
associated with eigenvalues near c. The procedure is as follows; 

(1) Use LCM to compute K - cM. 

(2) Use INV to factor K - cM. 

(3) Use EIG to compute the eigenvalues of eq. (9), (r - c), 
near zero. 

It is noted that INV can be used to determine the number of 
eigenvalues present in any range of interest. INV prints the number 
of negative terms encountered in factoring K - cM, which is equal to the 
number of roots below c in eq. (1). 

EIG measures convergence of the overall process by tracking 
changes in the eigenvalues of eq. (4), as computed in successive steps. 

Through RESET controls, the analyst may choose eigenvector initial 
approximations comprised of: 

(1) random vectors generated internally by EIG, and/or 

(2) data sets, from any source, resident in a library. 

If a vibrational eigenproblem is solved, output from EIG consists 
of two data sets: (1) VIBR U, the final approximations (whether converged 
or not) of all system eigenvectors being iterated on, and (2) VIBR EVAL,, 
the associated eigenvalue approximations ( not in HZ). If a buckling 
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eigenproblem is solved, the output is (1) vectors in STAB U, and 

eigenvalues in STAB EVA I,. Eigenvectors are normalized to X*MX = 1 

or X*K_ X = 1. 

G 


RESET Controls. 


Default 


K 


KG 


M 


MTYP 


CON 


K 


KG 


M 


SPAR 


Name 

Value 

Meaning 

PROB 

VIBR 

VIBR indicates solution of eq. (1). 
STAB indicates solution of eq. (2). 

NDYN 

5 

Maximum number of passes through the 
EIG iterative process. 

CONV 

. 0001 

See below. 

LOWER 

1 

See below. 

UPPER 

.8m 

For iteration j, the measure of convergence 
of eigenvalue i is e. = 


<r i j - r i j ' *> / *2 


EIG terminates when ej is less than 
CONV for eigenvalues LOWER through 
UPPER. The default value of UPPER is 
80% of the number of functions being used 
in the iterative process. 

First word of the name of the data set to be 
used as K in eq. (1) or (2). Common 
RESETS are K = K+KG, and K = K+CM. 

First word of the name of the data set to 
be used as Kg in eq. (2). 

First word of the name of the data set to be 
used as M in eq. (1), A common RESET is 
M = M+RM. 

Second word of the data set to be used as M 
in eq. (1). If element masses are to be 
ignored, RESET M = RMASS, MTYP = DIAG, 

Constraint case applied in factoring of K, or 
the K surrogate, by INV. 
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RESET Controls. 


Name 


Default 

Value 


KLIB 1 


KILIB 1 


KG LIB 1 

MLIB 1 

OUTLIB 1 


INLIB I 


Li 


Meaning 

Library containing K, or the K surrogate. 

Library containing the factored K (or K 
surrogate) produced by INV. 

Library containing KG, or the KG surrogate. 

Library containing M, or the M surrogate. 

Library into which EIG output is to be 
delivered. 

Library containing an existing data set 
(named either VIBR U or STAB U) from 
which vectors LI through L2 (see below) 
are to be extracted for use in the iterative 
procedure. 

See below. * 


L2 


See below. 


I NIT 0 


XRAND 0 

SHIFT 0 


In addition to, or instead of, initial vectors 
obtained from data sets in INLIB, EIG 
generates INIT vectors containing random 
numbers. That is, the total number of 
functions being iterated is generally m = 

INIT + L2 - LI + 1. If INIT = 0, and LI and 
L2 are not specified via RESET, EIG will 
use as initial vectors all vectors from VIBR U 
(or STAB U, if PROB = STABILITY), in INLIB. 
Note that it is allowed to have OUTLIB = INLIB, in 
which case results from previous executions are 
not^ overwritten; they are only disabled by the 
creation of a new data set of the same name 
(e. g. VIBR U), and can be re-enabled via DCU. 

Advance the random number generator by 
XRAND numbers (to avoid re-creating the 
sarne vector as in a previous execution). 

SHIFT is added to the computed eigenvalues 
(i.e. shift corresponds to c in eq. (9) ), 
before they are stored in VIBR EVAL or 
STAB EVAL. 
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RESET Controls. 


Default 

Name Value 

MAXIT 1 

ROUNDOFF .0001 


REPRT 0 

HIST 4 


TIME 0 

RRCH 0 


Meaning 

Same as in processor DSOE. 

Same as in processor DSOL, If 
RESET, ROUNDOFF and. CONV should 
be equal. 


Non-zero value results in printout of 
roundoff error measures during linear 
equation solutions. 

Non-zero values, 4 or 8, result in 
printout of eigenvalue convergence 
history, displayed as HIST columns 
per page. 

Non-zero values result in printout of 
intermediate CP and wall clock times. 

Non-zero values result in execution of 
a check of the accuracy of each Rayleigh 
Ritz analysis. 6 


Following any RESET cards, two additional types of commands may appeal 
(1) RSCALE Sj, s ^ , s 3 , s 4> Sg, s^. 

The above causes all direction - i joint motion components 
of all random vectors produced by the INIT RESET control 
to be multiplied by s.. Default values are s. = 1. for i = 1, . 

2, 3, and s. = . 01 for i = 4, 5, 6 (rotations). 


(2) PRINT K1 K2 K3 K4 K5 

Nonzero values of the K's result in production of the 
following categories of printout: 
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Kl: Vectors Y, KY , and MY are printed during each 

iteration. Bring a boxcar to haul your printout 
home . 

K2: Terms in the m x m coefficient matrices of eq. (4) 

are printed during each iteration. 

K3: All Rayleigh -Ritz eigenvectors, Q. in eq. (4), are 

printed during each iteration. 1 

K4: Vectors of eq. (5) are printed during each iteration. 

K5: The final set of eigenvector approximations are printed. 

Kl through K4 should be used only to diagnose abnormal result 

Suggested Technique . The RESET controls allow EIG to be used in 
many ways. The best approach for a particular case will depend on many 
factors, such as problem size, the number of eigensolutions required, the 
mode of execution (batch or interactive), etc. Accordingly no fixed rules 
of operation can be given; however, a few guidelines of fairly general 
applicability can be stated; 

(1) If eigenvectors for a similar structure are known, use them 
as initial approximations. This is very useful in parametric 
studies of the effects of changes in section properties, mass 
distribution, etc. 

(2) To compute modes of structures for which rigid body motion is 
possible, use the spectral shift method. The shift constant, 

c in eq. (9), should be selected within the expected range of 
eigenvalues associated with the elastic modes. Too small a 
value of c results in excessive roundoff errors due to 
K - cM being nearly singular. 
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(3) In the first execution of EIG, choose a small value for 

NDYN , perhaps as small as 2 or 3 for large structures. 
This will give the analyst useful information about the 
eigenvalue distribution in many cases, providing a basis 
for determining whether vectors should be added or deleted 
in the ensuing execution. For example, suppose you want 
to compute vibrational eigenvalues from . 0 to 100. ; and 
that you RESET INIT =12, NDYN =3. If it is found that 
only the lowest three eigenvalue approximations are less 
than 100. , it would likely be reasonable in the next run to 
continue iterating on only the lowest five or six modes; e.g. 
RESET LI = 1, L2 = 6. 

Generally, if p modes are to be computed, at least 1.2 p 
vectors should be used in the EIG iterative process, since 
the higher eigenvalues may converge much more slowly thar 
the mtenor (close to zero) eigenvalues, depending upon the 
spread of neighboring eigenvalues. 

If good turnaround is available (e.g. via an interactive 
terminal), it is almost always best to proceed in a sequence 
of short executions to avoid wasting computer resources by 
iterating on either too many or too few vectors. Note that 
the default operation ( @ XQT EIG, followed by no input cards 
is to resume iteration on all vectors computed and stored in 
a library in the last execution. 
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{4) If a very large number of modes is required, it may be 

cost effective to make several shifts and compute clusters 
of modes spanning the desired eigenvalue range. The 
resulting several data sets can subsequently be merged, 
if necessary. 

Core Requirements. In the following discussion, m is the number 
of vectors being iterated, L is the block length of the factored K (or 
K surrogate), and J is the number of joints. There are too parallel 
requirements, the second of which will dominate in almost all cases: 

(1) 3 m^ + I6n, and 

(2) 2 m + L + B, 

where B is a pool area usable only in integral multiples of 12J„ The 
10 activity of EIG will diminish linearly as B is increased in multiples 
of 12 J. 

Code Release Data . Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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Section 4 

DCU, THE DATA COMPLEX UTILITY PROGRAM 

Function - The term "data complex" refers collectively to all SPAR 
direct access libraries involved in an analysis. In defining DCU command 
syntax, the symbol Lib will be used repeatedly to represent a SPAR library 
internal designation (1, 2, etc. ), and Id will identify one or more data sets. 
Id may be in any of the three forms indicated below, unless specifically 
stated otherwise: 

(1) The four-word data set name, MASK-filled. 

Example: K SPAR is the same as 

K SPAR MASK MASK. 

(2) An integer, n, indicating the data set associated with 
sequence number n in the table of contents. 

(3) Integers n, m, indicating (m-n + 1) data sets associated 
with sequence numbers n, n+1, - - - m. 

(4) Omitted - meaning all data sets in a library. 

Consider as an example, the most commonly executed command, 

TOC Lib Id. 

Examples: 

• TOC 1 K SPAR$ Id form (1) causes single line of the 
Table of Contents of library 1 to be printed, i.e. the line 
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corresponding to the first data set named K SPAR MASK 
MASK to be located. 

• TOC 1 27$ Id form (2) causes printout of line 27 of the 

Table of Contents of library 1 . 

© TOC 1 32, 50$ Id form (3) causes printout of lines 32 through 

50 of the Table of Contents of library 1. 

® TOC 1 $ Id form (4) causes the entire Table of Contents of 
library 1 to be printed. 
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Other commands currertly available in DCU are summarized below. 

® ^ib ld$ . Data set(s) are marked as disabled. 

The data set(s) are still present in the library, but cannot 
be accessed until they are re-enabled via the following command: 

• gN/iBLE Lib, ld$ . Only forms (2) or (3) of Id are allowed 
for this command. 

• PRINT Lib, Id, i^, i j^, v^, v^. The identified data set(s) ar< 

prmted. Not all SPAR-generated data sets can be printed via 
this command. Only data such as tables, stresses, eigenvectors, 
displacement vectors, etc. , in regular "table" form should be 
printed in this way. Data such as topological maps (from TOPO), 
sparse matrices (from K, M, KG, INV, LCM, etc.), and element 
information packets (from E, EKS) should not be named on a 
PRINT command. 

The items, if given, must be the 7th, 8th, - - - items 

on the card. For eigenvector printout these parameters indicate component 
*1 through i 2 , joints jj through j 2# vector (blocks) Vj through v^ These iten 
have similar meaning for other multi-block data sets. 

• COPY Lib^ Lib 2> Id$. Copy the indicated data sets from Libj 
to Lib 2 . Disabled data sets are not copied. 


4-3 



DCU 


• XCOPY Libj, Lib 2> Id$. The indicated data (in Lib^) is 
written as an ordinary sequential file, Lib 2 , comprised of 
physical records identical to individual blocks of the data 
set as it resides in Lib^. As an example, XCOPY 1, 5, 
VIBR U$ causes eigenvectors to be written onto a file known 
externally as SPAR-E (UNIVAC) or SPARLE (CDC), cor- 
responding to Lib 2 =5. The output file will contain one 
physical record for each eigenvector. 


® SCALE C^, C 2 $. This command may precede a COPY 
or XCOPY command. If it does, the output data sets 
will = Cj (input) + C 2 . 


• CHANGE Lib Id o j d> Id new> This command causes the 

name of data set to be changed from Id .. to Id . Only 

old new y 


the full 4- word name form is permitted for both Id's. 


® DUPLICATE Lib 1 , Lib 2 $. Lib 2 is created identical to 
existing Lib^, including disabled data sets, if any. 


® TWRITE Lib$ . Lib is written onto tape nt .{see NTAPE 
command). The complete library is written in physical 
records as large as the available core will allow. 


® TREAD Lib$ . Lib is read from tape nt (see NTAPE command). 
Available working core space must be as large as it was when 
the TWRITE was executed. 
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® NTAPE = nt$ . The internal unit number of the tape to be 
used in the next TWRITE or TREAD command is nt 
(default=20). Note that logical 20 is Univac file SPAR-T, 
CDC file SPARLT. 

® STOR E Lib> ID$. Lib is stored as a data set named Id, 

in library nl (see LIB LIB command). Id may only be a full 
4 word name, the first two words of which are typeless. 

• RETRIE VE Lib, Id$ . The data set Id is recovered from the 
library nZ and constituted as library Lib. 

• LIB LIB = nX$. The internal unit number of the library 
library is n£ (default=12). Note that logical 12 is Univac 
file SPAR-L, CDC file SPARLL. 

® TITLE LIB 1 - - - Alphanumeric title for Lib - . 

The label -field title is imbedded in Lib, and will be 
displayed at the beginning of each table of contents 
printout produced by a TOC command. 

• STATUS Lib$. The number of library entries and the 
current I/O counts for Lib are printed. 

• ABORT n$. To cause an error-abort if an abnormal 
event occurs in DCU, set ABORT=l$. 
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Core Requirements . Working core must be sufficient to accommodate 
one block of each data set transmitted through core (e.g. via COPY, 
XCOPY, etc. ). See also the discussion of TWRITE and TREAD. 

Code Release Data . Level 7 (Univac, CDC), July 1974, coded by 
W. D. Whetstone. 
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Short List of Commands and Data Sequences 


TAB 


START 

UPDATE 

FORMAT 

NREF 


number of joints, list of excluded joint 
motion components. 

1 to enter update mode, 0 to leave. 

i, to select alternate format i. 

n, to identify reference frame n. 


MOD = m, to modify table entry numbers or 

joint numbers by m. 

After each sub-processor execution, 

TAB returns automatically to FORMAT=l, 
NREF-1, and MOD=0. 
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TAB Sub-processors 
TEXT: 


MATC: 

NSW: 

ALTREF: 

JLiOC: 


JREF: 

MREF: 

BRL: 


4/8 - - - 76 characters of text. 

k, E, Nu, Rho, Alpha, Stress. 

k, weight/ (area or length). 

k, ipa-p i 2 ,a 2 , i 3 ,a 3 , x 1 ,x 2 ,x 3 

k, xaj.xa 2 .xa 3 , xbj , xb 2 , xb 3 , ni, ijump, nj 

jjump, xcj, xc 2 , XC 3 , xdj , xd 2 ,xd 3 . 

Formats 1,2 = rectangular, cylindrical. 
NREF = reference frame. 


NREF= n: jj,j 2 ,inc: kj,k 2 ,inc: 

k. nb, ng, isign, c (FORMAT 1 ) 
k, ij, xy , x 2 , x 3 (FORMAT 2 ) 


NREF - n allowed under Format 2. 


k, 


.1 
i , 


1 

x l> 





x 2 
3 * 


A -2 



B A = BOX, k, bj.tj, b 2 ,t 2> 

TEE, k, b 2 ,t 2 - 

ANG, k bj,tj, b 2 ,t 2 . 

WFL, k, b 1 ,t 1 , b 2 ,t 2 , b 3 ,t 3 . 
CHN, k, b^tj, b 2> t 2> b 3 ,t 3 

ZEE, k, b^, t^, b 2 ,t 2 , b 3 ,t 3 . 
TUBE, k, inner r, outer r. 


GIVN, 

k, I x 

> o , 2> 

l2> &2 

1 > 3- > f 

t 

DSY, k 

' I I* 

oci, I 

z> <*Z> 

a, f, . 


V 

q 2> q : 

V y ll 

’ y 12’ " 

k ' S ll. 





S 21 ’ 

S 22* 




S 31 ’ 

s 32 J 

S 33. 



s 4r 

s 42' 

s 43’ 

S44. 


S 51 ' 

s 52’ 

s 53’ 

3 54’ 

s 55* 

s 6l » 

s 62> 

s 63> 

s 64> 

s 65» s 66 


BC: 

BD: 


^ , S , , S , 3, , a . « 

t c t c max mxn 
k, A, Ij, Qfj, h, c, q 1 


SA: 

k, t, t^, t 22 , t 12 . 

(Format= 


k ’ t ll’ t 22’ t l2’ 

(Format= 


c ll* c 12* c 22’ c 3 1 ’ 

C 32’ c 33* 


d ll’ d 12* d 22’ d 3 1’ 

d 32* d 33' 

SB: 

k, t, s, t , t 

max min. 



L* z i* z 2' 6> 

( card 1 ) 

' ‘ y 41 J y 42. (card 2) 


isotropic) 
aeolotropic) 
(membrane stiffness) 
(bending stiffness) 


A-3 



CON n: 


(constraint case n) 


JSEQ: 


RMASS: 


ZERO j ^ , j., 1 , inc*: j^ 2 , j^, 2 , inc 2 

NONZERO (same as above) 

RELEASE (same as above) 

SYMMETRY PLANE = n 
ANTISYMMETRY PLANE = n 
FIXED PLANE = n 
XZERO = x. 

RZERO — r, 

REPEAT n\ inc\ jmod^: 

^ 3 * J4/ J5 “ “ : 

■* 6 ' h’ ^ 8^9 : 

REPEAT n 2 , inc 2 , jmod 2 ; 

j 10’ j Il’ j 12^J 13 : 

etc. 

k, M, I r I 2 , i 3 . (Format 1 ) 

k, M^ , M^, M^ , Ij, I 2 * • (Format 2) 

REPEAT n, jump. 

CM = (M multiplier), ( I multiplier). 
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Sub -process or 

Name Element type 


E21 

E22 

E23 

E24 

E25 


General beam 

Finite length, directly specified K 
Axial element 
Plane beam 

Zero -length, directly specified K 


E31 

E32 

E33 


Triangular membrane 

Triangular bending element 

Triangular membrane + bending element 


E41 

E42 

E43 

E44 


Quadrilateral membrane 
Quadrilateral bending element 
Quadrilateral membrane + bending element 
Quadrilateral shear panel 


Pointer Default 
Name Value 

NMAT 1 

NSECT 1 

NOFF 0 

NNSW 0 


Associated Table 
Material constants (MATC) 

Section properties (BA - — — SB) 
Rigid link offsets (BRL) 
Nonstructural weight (NSW) 


NREF 1 


Beam orientation (MREF), Also 
indicates pressure sense for 3 and 4 
node elements. 
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Data Modifiers: 


MOD JOINT NUMBERS n, or MOD JOINT=n 
MOD GROUP NUMBERS n, or MOD GROUP=n 
MOD NMAT {or NSECT, etc.)=n 
INC NMAT (or NSECT, etc. )=n 

All MOD and INC parameters are reset to zero upon 
conclusion of each sub-processor execution. 

2 - node elements: 

•* 1 ’ ^2 

j ^ > j 2 » i* nj, jinc. (nj gridlines) 
jj> J 2 » ni, nj, jinc. {nj rings) 

j 1 * ^2’ ^ nj> jinc. (nj sets of spokes) 

3 - node elements: 

•V h’ ^3 

j]/ Jz’ J3’ n ^» n J' (rectangular mesh) 

j 1 ’ ^2' ^3* 11 *" n i» jinc. {nj 'fans'). 

4 - node elements: 

Jl> ^2’ ^3* J 4 

j 2 ’ ^3’ ^4' » n ^> kinc. (nk rectangular meshes) 

J 1 ’ * 2 * iy ^4’ 2 ’ ni> n J* {cylinder, nj bays) 
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Reset 

Name 

Default 

Meaning 


G 

1.0 

Mass = Weight/G 


IBEAM 

0 

If not zero, beam rotatory inertia ignored. 


INERT 

0 

If not zero, only displacement-dependent 
terms (as in bars, membranes) are 
included, for all element types. 

KG 





Reset 

Name 

Default 

Meaning 


QLIB 

i 

Library containing the prestress solution. 


QSEQ 

1 

Load case required number identifying 




the prestress solution. 

INV 

Reset 

Name 

D efault 

Meaning 


CON 

1 

Constraint case. 


K 

K 

Name of matrix to be factored. 


KILIB 

1 

Destination library for factored matrix. 


ILIB 

1 

Library containing AMAP (from TOPO). 

LrCM 

Forms R 

=C1 * Q1 + 

C2-*Q2 


Reset 

Name 

Default 

Meaning 


R 

K+KG 

Name of output matrix. 


Q1 

K 

Name of Q1 data set. 


Q2 

KG 

Name of Q2 data set. 
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Reset 

Name 

Default 

Meaning 

Cl 

1. o 


C2 

1. o 


RLIB 

1 

R destination library. 

Q1LIB 

1 

Q1 source library. 

Q2LIB 

1 

Q2 source library. 
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Q 


Reset 

Name 

Default 

Meaning 

QLIB 

i 

Destination library for Q output 

MLIB 

i 

Source library, system M 

M 

M 

System M, first word of name 

MTYPE 

SPAR 

System M, second word of name 

TMULT 

1 . 

T = TMULT * (T Input) - TNULL 

TNULL 

. 0 


PMULT 

1. 

P = PMULT * (p Input) - PNULL 

PNULL 

.0 



Commands : 

CASE n 1 - - - title for case n - - - 

NODAL FORCES: - - data - - 

NODAL MOTIONS: - - data - - 

NODAL TEMPERATURES: - - data - - 

NODAL PRESSURES: - - data - - 

INERTIA FORCE direction, acceleration (not g's) 

Nodal force, motion data cards : 

jo, k, q, ni, iinc, dqi, nj, jinc, dqj 


Nodal temperature, pressure data cards : 
jo, q, ni, iinc, dqi, nj, jinc, dqj 


DSOL 


Reset 

Name 

QLIB 


Default 


Meaning 

Source library for data from Q, 
destination library for output. 
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Reset 

Name 

Default 

Meaning 

CON 

1 

Constraint case 

LI 



L2 


Process load cases Li through L2 

K 

K 

Name of system K 

KLIB 

1 

Library containing system K 

ROUNDOFF 

10" 4 

Roundoff error tolerance 

MAXIT 

1 

Maximum iterative accuracy 
improvement passes 

EP 

1 

Nonzero: print error measures 

UP 

0 

Nonzero: print displacements 

FP 

0 

Nonzero: print reactions, nodal 
balance 


LCD 

Cases n, m' - - title for cases n through rn - - 
Lib^, seq^, increment^ 

Lib^, seq^, increment^ 

Lib^, seq^, increment^ 

GSF 

Reset controls QLIB, Li, L2 same as in DSOL 
Optional data cards: 

ELTS E23, E42, etc. 
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PSF 


Same as for GSF, Also reset LINES 
lines of printout/ page. 
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eig 


Reset 

Name 

Default 

Meaning 

PROB 

YIBR 

VIBR or STAB 

CON 

1 

Constraint case 

K 

K 


KG 

KG 


M 

M 


MTYPE 

SPAR 


ndyn 

5 

Max eigenvector iterations 

CONY 

io" 4 

Eigenvalue convergence test 

LOWER 

1 

Terminate when eigenvalues LOWER 

UPPER 

8 x no. fcts . 

through UPPER have converged 

LI 


Select vectors LI through L2 from 

L2 


previous solution for further 



iteration 

INIT 

0 

Generate INIT random vectors 

SHIFT 

. 0 

Add SHIFT to computed eigenvalues 

reprt 

0 

Nonzero: Print linear equation 
solution error measures 

MAXIT 

1 

Same as DSOL 

ROUNDOFF 

io" 4 

Same as DSOL 

Data Cards: 



PRINT 

k r k 2* k 3> 

k 4’ k 5’ ^ k 5 = f i n al eigenvectors) 

RSCALE 

s l’ S 2’ s 3’ 

S 4' S 5’ S 6 


Random vector element uj (direction j) is multiplied by s 

j 
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TITLE 

STATUS 

TOC 

DISABLE 

ENABLE 

COPY 

XCOPY 

SCALE 

CHANGE 

duplicate 

ntape 

T WRITE 
TREAD 

LIB LIB 

STORE 

RETRIEVE 

PRINT 

ABORT 


Lib Title - - 
Lib 

Lib, Id 
Lib, Id 
Lib, Id 


Lib^, Lib^, Id 
Libj, Lib 2 , Id 

Cj, C^ Output = Cj x (input) + C 2 


Lib ’ Id old- 
Libj , Lib^ 


Id 

new 


nt 

Lib 

Lib 


nl 

Lib, Id 
Lib, Id 

Lib, Id, i r i 2> j r j 2 , Vj , v 2 
n 
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ELEMENT FORMULATIONS 


Sections 1, 2, and 3 of this appendix summarize the assumptions involved in 
the routines used by SPAR to compute individual element K, M, and Kg matrices, 
respectively, relative to local element reference frames. The correspondence 
between these routines (named K21, TM, H3k, GQM, etc.) and the SPAR element 


type designations (E21, E22, etc.) is 

indicated in 

Table B-l. 



TABLE B-l SUMMARY OF 

K, M, and Kg 

ROUTINES 

Element 




^oe 

— 

K 

M 

M 

E21 

Gen. beam 

K21 

M 62 

BEAMKG 

E22 

Direct K 

(see text) 

none 

BARKG 

E23 

Bar 

(see text) 

M3 2 

BARKG 

E2li 

PI. beam 

(see text) 

M62 

BEAMKG 

E25 

Direct K, 

L=0 (see text) 

none 

none 

E31 

Membrane 

TM 

M3 3 

GTM 

E32 

Bending 

TPB7 

M63 

none 

E33 

Mem.+bend. 

TM+TPB7 

M63 

GTP 


Elil 

Membrane 

QMB5 

M3ii 

GQM 

E^2 

Bending 

QPB11 

M64 

none 

EU3 

Mem. +bend. 

QMB3+QPB11 

m61i 

GQP 

El \h 

Shear Panel 

QMB1 

M3ii 

none 
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1. Stiffness Matrices 


Individual element K's computed by the routines K21, TM, QPB11, etc. indicated 
in the first column of Table B-l are relative to intrinsic reference frames 
imbedded in the elements. As indicated on Figure B-l, the intrinsic frames 
— ~y e the_ ele ments as the structure deforms . Some of the principal reasons 

for using intrinsic frames are summarized below. 

(1) K's relative to intrinsic frames totally represent each element's inherent 
flexibility characteristics in minimal (least data) form. 

(2) Transformation, assembly into system K, etc. , can be performed by very 
general routines which are totally unaffected by the details of specific 
formulations used in computing the element intrinsic K's. As a result, 
the cost and time required to develop routines to implement various 
element formulations are minimized, since their total function is to compute 
an intrinsic K, given the dimensions, section properties, and material 
constants of the element. 



(1) Dashed lines indicate undeformed shape. 

(2) For beams, the origin of the intrinsic frame is imbedded in (rotates with) 
joint J2 > to which the beam terminus connects. 

(3) For bending elements, the origin of the intrinsic frame is imbedded in Jl. 
(ii) In the undeformed state, the intrinsic frames coincide with the local 

"element reference frames." 

FIGURE B-l INTRINSIC FRAMES 
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For each class of element, an intrinsic deformation vecto r, U, is defined as 
follows, 


Beams ; 

U 

3 -node 


membranes : 

U 

ii-node 


membrane s : 

u 

3 -node 


bending ; 

u 

ii-node 


bending : 

u 

In the above, the 

u J . 


11 1 1 1 1 

U 1 U 2 U 3 9 1 e 2 9 3 

2 3 3 it 

U 1 U i U 2, 

2 3 3 k U 1 t 

U 1 U 1 U 2 U 1 U 2 i 


' Kj ® 2 1 9 \ 9 2 I' 




‘i 


(i-D 


P * j (i.e. Jl, J2, . ..) in direction i, relative to the intrinsic frame 
For each type of element, the associated intrinsic stiffness matrix, K, is 
defined by the following equation , 


of 


Strain energy = V = 1/2 U t K U . 


( 1 - 2 ) 


The order of K is 
indicated above. 


6, 3, 5, 6, and 9, respectively, for the five classes 
For linear elastic materials, K is symmetric. 


of U 


The methods used to transform and assemble 
etc., will not be discussed here, since the 
to explain the basis of the intrinsic K's. 


the intrinsic K's into system ma trice 
purpose of the present discussion is 
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1«1 Beam Intrinsic Stiffness Ite trices 


For element types E22 and E25, the 6x6 intrinsic K, as defined by equations 
(1-1) and (1-2), is directly read as input. The following discussion defines 
the basis of E21 and (as special cases), E23 and E2l*. The characterizing 
feature of the E21 class of elements is that the force-deflection behavior may 
be expressed in the form indicated by equation (1-3). The relation between 
the <J«s and the intrinsic deformation vector, U, will be established subsequently 


A' 


S n 0 

0 

0 

S 15 

0 


H 
L 

CM 


s 22 

0 

S 2k 

0 

0 


4 2 1 

f 3 



S 33 

0 

0 

0 


*3 1 

f u 




S kh 

0 

0 

: ^ 

f s 

f 6 


symmetric 



8 SS 

0 

s 66 


56 


(1-3) 


( 1.-10 


In the above, the f ± 's are forces (i=l,2, 3) and moments (i=U,£,6) acting on a 
PiS-ldplan# imbedded in the beam origin, with the beam terminus frbced. The 6 . 's 
are corresponding displacements and rotations. As defined on Figure B-2, the 
transverse forces and bending moments do not act through the beam origin; instead 
they act in principal bending planes - hence the form of equation (1-3). 


* The rigid plane is perpendicular to the 3-axis of the intrinsic frame. 



The 1, 2 and 3 axes are those 2 
of the element, intrinsic 
frame . The l' and 2' axes 
are in cross-section principal 
planes . 




0 = angle of inclination 
of section principal 
axes relative to 
element intrinsic 
frame . 



Centroid 


The <5's are displacement 
and rotation components 
corresponding to the f ! s. 


FIGURE B-2 


FORCE AND MOMENT RESULTANTS EXERTED ON THE 
MEMBER END SURFACE BY THE TRANSVERSE RIGID PLANE 



In equation (1-3), s.^, snd are associated with bending about principal 
axis 2; s^, s^, and with bending about principal axis-1. 

On Figure B-3, components of the element intrinsic deformation vector, and the 
associated point forces, p , and moments are shown. 



In order to identify K, the equilibrium relation between the p . 's (Fig. B-3) 
and f^s (Fig. B-2) will be established'. In the following, 9, and z are as 

defined on Figure B-2. Also, c = cos 9, s = sin Q. 


K 


p 4 

! p 5! 

W 


or. 


1 

0 

0 

0 

0 

l-z 


0 

1 

0 

0 

0 

z_ 


0 

0 

1 

0 

0 

0 


0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

1 

0 


0 

0 

0 

0 

0 

1 


c 

s 

0 

0 

0 

0 


-s 

c 

0 

0 

0 

0 


0 

0 

1 

0 

0 

0 


0 

0 

0 

c 

s 

0 


P « T f, 

where T is defined by equation (1-5), and p - ^ p 2 . . . p^j 1 . 


0 

0 

0 


-s 


c 

0 


0 

0 

0 

0 

0 

1 


( 1 - 6 ) 


The kinematical relation indicated by equation (1-7) may readily be verified. 

(1-7) 


.(1-5) 


<J = T t U 
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From equations (1-U), (1-6), and (1-7), 
p = T S T^U, or 

K = T S l t . (1-8) 

If the shear center and centroid coincide, and the intrinsic frame axes are in 
principal planes, 

K = S. 


The following symbols will be used in identifying the eight non -zero s. . »s. 


L 

A 

E 

G 



length of beam, 
cross sectional area, 
modulus of elasticity, 
shear modulus, 

cross-section moments of inertia about principal axes 1 and 2, 
respectively, 

transverse shear deflection constants (see Timoshenko ''Strength 
of Materials," Part 1, p. 170), 
uniform torsion constant, and 

non-uniform torsion constant (see Timoshenko "Strength of 
Materials," Part. 2, p. 255-273). 


The axial spring constant is 
s ^ = AE/L. 

The torsional constant is 



which assumes no warping of end cross sections . 

Under the usual assumptions of Timoshenko beam theory, 


<5. = f 


1 3 El, 


- f 


5 2 El. 


« -f 


1 2 El, 


+ f. 


ei 2 5 


«1 L 

GA 
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or where e n 


EE, 


-s "i andk i 5 


El, 


ei + LV12 ’ 


V 


(±_ + 6 l\ 1 


f. 



\3 T) ‘T 


i 


L 2 




1 

E i 

1 1 


4 > 


1 

L ^ ^ 


f 5 

L J 


and 


$ 


1 _ 

L 


_1_ 

2 


+ (t*?) 


From the above, the first set of bending stiffness constants may be identified as 


Also, 


s n 

■ E k i 

s l£ 

1 1 

HjcvJ 

It 

s 55 

/ l 1 
■(r- 

S 22 

= i k 2 

S 2U 

- -h: 


■ G* 


) k 2’ 


where the definitions of and e^ and analogous to those of and e^. 
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1®2 Membrane and Bending Elements 


In the displacement formulation of linear finite element networks, the boundary- 
motion of each element is written in terms of certain generalized displacement 
quantities , us u al ly the displacement and rotation components of system Joints 
to which the element is connected. Elements are said to be conforming if the 
relations between the system generalized displacements and element edge motions 
are such that displacements are continuous across element boundaries. The 
basic problem in determining the stiffness matrix of a conforming element may 
be stated as follows : given the boundary motions of the element, determine the 
interior stress state. Since the late 1950's, this problem has most commonly 
been approached through the assumed displacement field - minimum potential 
energy method. In the July, 1961; issue of the AIM Journal, T.H.H. Pian 
proposed a simple, general, and highly effective alternative: the assumed 

stress field - minimum complementary energy method, which is outlined below. 


In terms of generalized displacements, q, (e.g. joint motions), the strain 
energy of an element is 


V = h q K q. 

The minimum complementary energy principal may be stated as 
n c - ^ ~ J Uj_ ^ d(Area ) = minimum 


( 1 . 2 - 1 ) 

( 1 . 2 - 2 ) 


Area 


where V is strain energy in terms of stress components, a . u and S are 

xj x l 

direction i components of prescribed boundary motion and boundary force on the 
exterior surface. 


To apply the principal, the stress distribution, a, is written in terms of m 
undetermined stress coefficients, b h , . . b , as follows. 

J- ci W 


a 

B 

a 


P B, where 

V V" b m 

* 11 ’ ° 2.2 •** 


(1.2-3) 


and elements of P are functions of position coordinates, selected such that 
equations of equilibrium are identically satisfied within the element, regardless 

of the values of the b . * s . 

x 
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Where the stress-strain relation is 
€ = N o , 

the internal strain energy within the volume 

V = j a 1 N <7 d(Vol) 

Vol 



where elements of L are functions of position coordinates, determined by the 
location of the boundary. 

Where S = [ s-^ s^ the relation between surface forces and the stress 

field coefficients is of the form 


S = E E, 

where the elements of R are functions of position coordinates determined by 
the location of the boundary, and by the form of P. 


Accordingly, 

* — r 1 _ t 

n c = — B HB-BTq, 

where 

T = J' R^ L d(Area). 

Area 

For minimum 77 , i.e. dfJ/db. » 0, for i = 1,2,..., 

C Cl 

HB = T q, 

B = H " 1 T q. 


(1.2-7 ) 


( 1 . 2 - 8 ) 
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and, from equation (1.2-U), 

V = h B* H _1 T B, 
from which it is concluded that 
K = H _1 T. 


(1.2-9) 


The procedure outlined above was proposed by Pian as a method of calculating 
K's approximating those of conforming elements. Where K is the stiffness 
matrix corresponding to the perfect conforming element, K approaches K 
from the soft direction (too much strain energy), as the number of assumed 
stress functions is increased. In the assumed displacement field — minimum 
potential energy approach, K approaches R from the stiff side as the number 
of assumed displacement functions is increased. 

A network comprised of conforming elements will converge uniformly from the 
stiff side toward the exact solution, as the mesh is made progressively finer. 
Even if a mesh were comprised of perfect conforming elements (equilibrium, 
compatibility equations satisfied over the interior of each element), the net- 
work solution would be too stiff. Therefore, K's based on the assumed stress 
field procedure may be viewed as a means of systematically introducing compen- 
sating error (K's softer than K), to obtain more accurate solutions. 

To summarize, the following are the two basic assumptions involved in applica- 
tion of the assumed stress field - minimum complementary energy method. 

(1) The boundary motion of the element is assumed as a function of system 
generalized displacements (e.g. joint displacements, rotations, etc.), 
in a form which provides inter-element compatibility. This amounts to 
defining the matrix L in equation (1.2-6). 

(2) An equilibrated stress field, expressed in terms of m undetermined coef- 
ficients, b-pb 2 , ... b > is assumed. This amounts to defining the 
matrix P in equation (1.2-3). 

To implement the above assumptions, the principal computational tasks are the 
calculation of H from equation (1.2-5), and T from equation (1.2-7). From H 
and T the element stress and stiffness matrices, H _1 T and T^H S are readily 
determined. 

The in-plane edge motion assumed for the TM (triangular membrane), QMB5 
(quadrilateral membrane), and QMB1 (quadrilateral shear panel) elements is 
shown on Figure Bl.2-1. The edge motion assumed for the TPB7 (triangle), 
and QPB11 (qua dr ila t era 1 ) bending elements is shown on Figure Bl.2-2. 
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L = length of edge 


deformed edge 
(straight line) 



V 4 k ♦ fark - 1 ’ 2 


FIGURE EL. 2-1 IM-PLANE EDGE MOTION ASSUMED IN TM, QHB5, QMB1 



w is normal to the element surface . 
w = a cubic function of s such that 


At node i. 


w = w 


<?w 

as 





At node 
w = 


3 > 

aw _ 
as 


z = normal to 



a t edge ij s 

<f> = 4* 1 + -|- ) 


FIGURE Bl.2-2 EDGE MOTION ASSUMED IN TPB7, QjPBll 
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Stress fields assumed for TM, QMB1, and QMB5 are indicated below. In the 

following, x and y correspond to directions 1 and 2 of the intrinsic frame 

defined on Figure B-l, and N , N and N are in-plane stress resultants. 

j «*y 

For TM (constant-stress triangular membrane). 


N = b 
x 3 


N = b 
y 2’ 


V 


= b 


y 


For QMB1 (quadrilateral shear panel). 


N = 0, 

x * 


N 


y " 

V ‘ b r 

Far QMB^ (quadrilateral membrane), 

"x ' b i * V- 

\ = b 3 + V' 

N “ V 

xy 5 

In TPB7 and QPB11, the (Kirchoff) stress resultant -curvature relation is of 
the following form. 


M 


M 

y 

M 

xy 


D 11 

*12 

D 13 

D 12 

D 22 

D 23 

D 13 

D 23 

D 33 


j 


w/dx 


2 , 2 
■d vt/dy 


2 d w/dx&y 


(1.2-10) 


The stress fields assumed in TPB7 and QPB11 are derived, using equation (1.2-10), 
from the displacement functions indicated below, which satisfy both equilibrium 
and compatibility* equations. 


This is not essential to application of the assumed stress field— minimum 
complementary energy method. It is necessary only that the assumed stresses 
satisfy all equations of equilibrium. 
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For TPB7 (triangle): 


w = “ b i“5“ “ b 2 TT + b 3 ^ ~ b ‘ ^ * 

For QPB11 (quadrilateral): 


\ “ 5 " 


'5 ^ " b 6 “2 


** v “ b 7 ^ (lo2-ll) 


w = same as for TPB7 

o 

H 

& 

[4 * °3 

where 


a l = 

- (2D 12 + 1 

°2 

2D 23 /D H j ’ 

a 3 = 

2 V d ii> 

II 

“ D 22 /D ll* 


8 1 


v 2 2 

X_j_ 


+ a -dL_ 

2 1 12 


xp 

J 9 \~T 


y \ 

h IT r 


jt 

2 12 


( 1 . 2 - 12 ) 


Calculation of H and T (see equations (1.2-5) and (1.2-7) for the above 
elements involves area and line integrals, respectively, of polynominals 
of varying degrees of complexity. For the simple terms, closed form expres- 
sions are readily determined. For the complicated forms, Gaussian integra- 
tion of order sufficient to assure exact integrals is used. 

The routines which compute intrinsic K's for the TM, QMB1, QMB5, TPB7, and 
QPB11 elements were written by C. L. Yen. 



2. Element Mass Ife trices 


For immediate reference, the SPAR element reference frame convention is shown 
below . 


2 


1 



2 -node elements 


J l> J 2 > J-jj all are in plane 1-2 
of the element reference frame 



FIGURE B2-1 ELEMENT REFERENCE FRAMES 


In the following, u* k and 0^ are the total displacement and rotation, res- 
pectively, of node J ± in direction k, relative to the element reference frames. 
Local generalized displacement vectors, q, for individual element categories 
are defined as follows. 

Where u 1 = (u 3 ^ u 3 ^ u 1 ^, 

for bars (M32), 

f 1 2l t 

q = [u u J . 

For triangular membranes (M33), 

[ 1 2 3 1 t 

q = u u u 


For beams (M52), 

q - [u\ u 1 ., U X 3 9^ e X 2 9 1 u\ u : 


2 * 2 _ 9 2 e 2 e 2 

2 3 1 2 3 


Where ^i = |u%, 0“ 0 


1 0 1 
1 y 2 


for triangular bending elements (MTP), 

Q - K 4 ? 4 ?}\ 
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For quadrilateral bending elements (MQP), 
q = [/ <A 2 ^ ‘ 

Where q is as defined above, the kinetic energy of an element is expressed in 

+■ 

terrns of a local element mass matrix, M, as h q Mq, Transformations to 
system coordinates, assembly into system M, etc., are performed in SPAR by 
general purpose routines which are independent of the source of the element 
M's. 


In the following discussion, x, y and z are position coordinates in directions 
1, 2, and 3 of the element reference frame; and u, v and w are displacement 
components in directions x, y, and z. In order to compute element M's, dis- 
placement functions in terms of n undetermined coefficients, a n a a 

1’ 2 s n* 

are assumed in the following form. 


n 


~v 

V 

= G 

a .z 

* 

« 

_ w - 


a 

L n J 


(2-1) 


where elements of G are functions of x, y, and z. Where m is mass/volume, 
the kinetic energy of the element is 


= h f m (u 2 + v 2 +w 2 )d(Vol). 


Vol 

Where A = |a-^ a g . . . a j t , substitution of (2-1) into (2-2) gives 


( 2 - 2 ) ' 


T = 


*1 / 
1 Vol 


m G G d(Vol) A 


(2-3) 


For each element type, n is selected equal to the number of elements in q, so 
that a unique kinematic relation between A and q of the form indicated by 
equation (2-10 is readily determined. 

q = GA > (2-U) 

where the elements of C depend on the locations of the node points. Substitu- 
tion of equation (2-10 into (2-3) identifies the element mass matrix as 


M = (C -1 ) 1 [ m G^ G d(Vol) 



(2-5) 
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The displacement fields assumed for each type of element are indicated below. 
For bars (M32), 


u = ^ + a 2 ,, 

V = a 3 + a k Z ’ 

w = + a^z. 


(2-6) 


For triangular membranes (M33), 
u = a 1 + a 2 x + a 3 y, 
v = 3^+3^+ a 5 y, 

w = a 7 + a 8 x + a ? y. 


(2-7) 


For beams (M62), where u, v, and w are displacements of the central 
axis, and where 9 is the total twist angle, 

U = a i + a 2 Z + a 3 z2 + > 

v - a s + a 6 s + a ?z 2 + a^ 3 , 

W = a 9 a lO Z> 

0 ' 3 11 * a i2 Z ‘ 


(2-8) 


For triangular bending elements (MPP) 

2 

w = + a 2 * + aye 

+ \y + + a 6 x 2 y 

2 2 2 2 
+a?y + a 8 xy + aye y . 

For quadrilateral bending elements (MQP) 

, _ „3 


a l 

+ a gX + 


+ agxy + 

2 

^io^ 2 

^ 3 

3 

^n y 

+ a 12 Xy * 


(2-9) 


( 2 - 10 ) 
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Using the displacement fields indicated above, the integration indicated by- 
equation (2-5) ) is readily performed in closed form, except for MTP and MJP, 
in which cases Gaussian integration is used. In the case of M62, the rotatory- 
inertia terms associated with all three cross section moments of inertia may- 
be included, if desired. 

The correspondence between the above and the designations used in Table B-l is 
as follows: 

M32: As defined above. 

M33 i As defined above . 

M31+: Two overlaid pairs of M33's. 

M62: As defined above. 

M 63 i M3 3 + MTP. 

M6U: M3U + MQP. 

It is noted that the above are generally not "consistent" mass matrices, 
according to the usual definitions since they are derived on the basis of 
arbitrarily selected displacement fieldswhich only in special cases coin- 
cide with the displacement fields associated with the corresponding element 
K's. They usually give a much better approximation of system kinetic energy 
than do lumped mass models, however. 

The routines which compute the local element mass matrices designated above 
as M33, M3U, M63 and M5U were written by C. L. Yen. 
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3. Geometric (Initial Stress) Stiffness Matrices 


The meaning of the initial stress stiffness matrix, K , may be viewed in the 

S 

following way. Assume that a structure in a given equilibrium configuration 
undergoes an additional deformation, q. The pre-stress dependent quadratic 


The 


terms in the associated change in strain energy are U = q K-q. 

it & 

ordinary quadratic terms are K q. Accordingly the total change in strain 

energy associated with the perturbation, q, is approximately h q t (K^+K)q. 

This provides means of (1) including effects of pre-stress on deformation of 

linear systems, and (2) performing piecewise linear analysis of geometrically 

nonlinear problems. The application of K to linear bifurcation buckling is 

based on the following: if K corresponds to an applied force system, F, the 

change in strain energy due to perturbation of a structure to which a load AF 

is applied is q^( XL +K)q = 0. For stationary strain energy, (\K +K)q = 0. 

G g 

In the following discussion of the basis of individual element K^'s, the 

notation (u,v,w, etc.) used will be the same as was used in the discussion 

of element mass matrices. 

For 2-node elements, the principal contributions to K are due to changes in 

s 

axial strain due to lateral deformation, u and v, 

\2 




Where P is the pre-stress axial force in the element, the corresponding 

z 

change in strain energy within the element is 
.L _L 


U-. 


- / 


P € d Z 
z z 


f 


/any 

\dzf 


m 


dz. 


(3-D 


A planar example is shown on Figure B3-1. 
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K g EFFECTS in a plane beam 


figure B3-1 



I $ 

Similarly , in two-dimensional elements the changes in membrane strains due to 
lateral displacement produce the following change in strain energy. 



+ n + 2N ~ 

y\dy xy dx dy 


d (Area ) 


(3-2) 


In the above, N , If and If are the pre-load membrane stress resultants . 


The procedures used to calculate are similar to those used in calculating 
element mass matrices. Lateral displacement functions are assumed in ternns 
of n undetermined coefficients, a^ a g , . . . a n . Except for the quadri- 
lateral membrane, GQM, each of the formulations assumes the same dis- 
placement field as a corresponding mass matrix, as tabulated below. 


K g lype 


Assumed 

Displacements 


BARKG 

BEAMKG 

GTM 

GTP 

GQP 


Bq (2-6), axial element. 

Bq (2-8), general beam. 

Bq (2-7), triangular membrane. 

Bq (2-9), triangular bending + membrane. 

Eq (2-10), quadrilateral bending + membrane. 


For the quadrilateral membrane, GQM, lateral displacement is assumed in the 
following form. '' 


W “ a l + a 2 X + a 3 y + 


The membrane stress distribution used in computing K for quadrilaterals GQM 
and GQP Is the same as indicated for the QMB5 element (i.e. linearly varying). 
BEAMKG includes, in addition to the principal terms indicated in eq. (3-1), 
other secondary terms involving initial shear stresses and bending moments. 

As before, the integrations indicated by equations (3-1) and (3-2) in some cases 
are readily determined in closed form; and in others are evaluated using Gaussian 
integration. These routines were coded by C. L. Yen. 


/ 

/ 
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